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Understanding the predictions of general relativity for the dynamical interactions of two black holes 
has been a long-standing unsolved problem in theoretical physics. Black-hole mergers are monu- 
mental astrophysical events, releasing tremendous amounts of energy in the form of gravitational 
radiation, and are key sources for both ground- and space-based gravitational-wave detectors. 
The black-hole merger dynamics and the resulting gravitational waveforms can only be calculated 
through numerical simulations of Einstein's equations of general relativity. For many years, nu- 
merical relativists attempting to model these mergers encountered a host of problems, causing 
their codes to crash after just a fraction of a binary orbit could be simulated. Recently, however, 
a series of dramatic advances in numerical relativity has allowed stable, robust black-hole merger 
simulations. This remarkable progress in the rapidly maturing field of numerical relativity, and 
the new understanding of black-hole binary dynamics that is emerging is chronicled. Important 
applications of these fundamental physics results to astrophysics, to gravitational- wave astronomy, 
and in other areas are also discussed. 
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I. PRELUDE 

The final merger of two black holes in a binary sys- 
tem releases more power than the combined light from 
all the stars in the visible Universe. This vast energy 
comes in the form of gravitational waves, which travel 
across the Universe at the speed of light, bearing the 
waveform signature of the merger. Today, ground-based 
gravitational-wave detectors stand poised to detect the 
mergers of stellar black-hole binaries, the corpses of mas- 
sive stars. In addition, planning is underway for a space- 
based detector that will observe the mergers of massive 
black holes, awesome behemoths at the centers of galax- 
ies, with masses of ~ (10 4 — 10 9 )Mq, where Mq is the 
mass of the Sun. Since template matching forms the ba- 
sis of most gravitational-wave data analysis, knowledge 
of the merger waveforms is crucial. 

Calculating these waveforms requires solving the full 
Einstein equations of general relativity on a computer in 
three spatial dimensions plus time. As you might imag- 
ine, this is a formidable task. In fact, numerical rela- 
tivists have attempted to solve this problem for many 
years, only to encounter a host of puzzling instabili- 
ties causing the computer codes to crash before they 
could compute any sizable portion of a binary orbit. 
Remarkably, in the past few years a series of dramatic 
breakthroughs has occurred in numerical relativity (NR), 
yielding robust and accurate simulations of black-hole 
mergers for the first time. 

In this article, we review these breakthroughs and 
the wealth of new knowledge about black-hole mergers 
that is emerging, highlighting key applications to astro- 
physics and gravitational-wave data analysis. We focus 
on comparable-mass black-hole binaries, with component 
mass ratios 1 < q < 10, where q = Mi/M 2 and M\, M% 
are the individual black-hole masses. We will frequently 
also refer to the symmetric mass ratio 
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M X M 2 



{Mi+M 2 f {l + qf 



(1) 



For simplicity, we choose to set c = 1 and G = 1; with 
this, we can scale the dynamics and waveforms for black- 
hole binaries with the total system mass M. In particu- 
lar, we can express both length and time scales in terms of 
the mass, giving M~5x 10~ 6 M/M Q s ~ 1.5M/M Q km. 



We begin by setting both the scientific and historical 
contexts. In Sec. [TT] we provide a brief overview of astro- 
physical black-hole binaries as sources for gravitational- 
wave detectors. We next turn to a historical overview 
in Sec. IIII1 surveying efforts to evolve black-hole merg- 
ers on computers, spanning more than four decades and 
culminating with the recent triumphs. Having thus set 
the stage, we focus on more in-depth discussions of the 
key components underlying successful black-hole merger 
simulations, discussing computational methodologies in 
Sec. IIV1 including numerical-relativity techniques and 
black-hole binary initial data. Section [V] is the heart of 
this review. Here we discuss the key results from numer- 
ical relativity simulations of black-hole mergers, follow- 
ing a historical development and concentrating on the 
merger dynamics and the resulting gravitational wave- 
forms. These results have opened up a variety of excit- 
ing applications in general relativity, gravitational waves, 
and astrophysics. We discuss synergistic interactions 
between numerical relativity and analytic approaches 
to modeling gravitational dynamics and waveforms in 
Sec. IVH and applications of the results to gravitational- 
wave data analysis in Scc. lVIll The impact of merger sim- 
ulations on astrophysics is presented in Sec. IVIII1 which 
includes discussions of recoiling black holes and potential 
electromagnetic signatures of the final merger. We con- 
clude with a look at the frontiers and future directions of 
this field in Sec. HX1 

Before we begin, we mention several other resources 
that m a y inte rest our readers. The review articl es by 
iLehnerl ( 200lh and iBaumgarte and Shapirol ( 20031 ) pro- 
vide interesting surveys of numerical relativity several 
years before the breakthro ughs in bla c k-hole merger sim- 
ulations. The article by iPretoriusI ( 20091 ) is an early 
review of the recent successes, co vering some of the 
same topics that we discuss here. iHannanJ (2009) re- 
viewed the status of black-hole simulations producing 
long waveforms (including at least ten cycles of the dom- 
inant gravitational-wave mode) and their application to 
gra vitational- wave data analysis . Fin ally, the tex t books 
by iBona and Palenzuelal (|2005l) and lAlcubierrd (|2008l ) 
provide many more details on the mathematical and com- 
putational aspects of numerical relativity than we can 
include here, and serve as useful supplements to our dis- 
cussions. 



II. BLACK-HOLE BINARIES AND GRAVITATIONAL 
WAVES 



Black holes and gravitational waves are surely among 
the most exotic and amazing predictions in all of physics. 
These two offspring of Einstein's general relativity are 
brought together in black-hole binaries, expected to be 
among the strongest emitters of gravitational radiation. 
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A. Basic Properties 

We begin by presenting some basic properties of black 
holes and gravit ational waves. For fu ller discussions an d 
more details, see lMisner et all (|l973l ) and lSchutd (|2009h . 



1. Black-Hole Basics 

A black hole forms when matter collapses to infinite 
density, producing a singularity of infinite curvature in 
the fabric of spacetime. Each black hole is surrounded 
by an event horizon, at which the escape velocity is the 
speed of light. The event horizon is a global property of 
the spacetime, since it is defined by the paths of "out- 
going" photons that are the boundary between photon 
trajectories that must fall inward, and those that can es- 
cape to infinity. The photons defining the event horizon 
hover at finite radius at the surface of the black hole. 
Since, in principle, mass (energy) can fall into the event 
horizon at late times - which will move the location at 
which photon paths can hover - we must know the en- 
tire future development of the system to locate the event 
horizon. 

When black holes merge, a single event horizon forms 
whose area is at least as large as the sum of the individ- 
ual horizons. Since numerical relativists want to know 
when this occurs during the course of a calculation, they 
rely on a related concept known as an apparent hori- 
zon, whose location depends only on the proper ties of 
the spacetime at any given time ( Poissonl 120041 ). For 
quiescent black holes, the apparent and event horizons 
coincide; for more general holes, the apparent horizon is 
always inside the event horizon (with restrictions on the 
behavior of the matter involved). So, in terms of causal- 
ity in a numerically-generated spacetime, any physical 
phenomenon found inside an apparent horizon should not 
leak out and affect the spacetime outside. 

The simplest black hole is nonrotating and is described 
by the spherically symmetric Schwarzschild solution to 
the Einstein equations of general relativity in vacuum 
(i.e., with no "matter" sources in the spacetime). A 
Schwarzschild black hole is fully specified by one quantity, 
its mass M. The horizon is located at coordinate r = 2M 
(in Schwarzschild coordinates); its area is 4ir(2M) 2 . 

More general black holes can have both charge and 
spin. Since a charged black hole in astrophysics will gen- 
erally be neutralized rapidly by any surrounding plasma, 
we can consider only rotating, uncharged black holes. 
Stationary (i.e., time independent) black holes are de- 
scribed by the axisymmctric Kerr solution. A Kerr black 
hole is fully specified by two quantities, its mass M 
and its angular momentum per unit mass a. The event 
horiz on is located at the Boyer-Lindquist ( Misner et all . 
11973b radius r+, where 



Equation © requires a < M\ when a = AI the black 
hole is said to be maximally rotating or "extremal" . No- 
tice that r+ = 2M when a = 0, and that r+ decreases 
as a increases, thus bringing the location of the horizon 
deeper into the potential well as the black-hole spin in- 
creases. 

Photons and test particles in the vicinity of a single 
black hole can experience either stable or unstable orbits. 
For a Schwarzschild black hole of mass M, the innermost 
stable circular orbit (ISCO) occurs at r = 6AI for massive 
test particles. In the case of a Kerr black hole, the ISCO 
is closer in for co-rotating test particles, and farther out 
for counter-rotating particles. 

While the concept of an ISCO is strictly defined only 
for massive test particles, it has proven useful for stud- 
ies of the spacetime around two black holes spiralling 
together on quasicircular orbits. Imagine that you put 
the two black holes on an instantaneously circular orbit 
around each other; at that moment they have neither 
nonzero radial velocity nor nonzero radial acceleration. 
At any given separation, the black holes have some an- 
gular momentum. The ISCO is the separation where 
that angular momentum is a minimum, in analogy to 
the test particle definition. Black holes at closer separa- 
tions would be expected to fall inward, toward the center, 
even without radiating angular momentum via gravita- 
tional radiation. 



2. Gravitational Wave Primer 

Gravitational waves arc ripples in the curvature of 
spacetime itself. They carry energy and momentum and 
travel at the speed of light, bearing the message of dis- 
turbances in the gravitational field. 

As with electromagnetic waves, gravitational waves 
can be decomposed into multipolar contributions that 
reflect the nature of the source that generates them. Re- 
call that electromagnetic radiation has no monopole con- 
tribution due to the conservation of total charge. By 
analogy, conservation of total mass-energy guarantees 
that there can be no monopole gravitational radiation. 
Since dipolar variations of charge and currents are pos- 
sible, electromagnetic waves can have a dipole character. 
However, conservation of linear and angular momenta re- 
moves any possibility of dipolar gravitational waves, so 
the leading-order contribution to gravitational radiation 
is quadrupolar. 

Gravitational waves are thus generated by sys- 
tems with time-varying mass quadrupole m oments 
( Flanagan and Hughesl feoOS; Mi sner et In the 

wave zone, a gravitational wave is described as a pertur- 
bation ft to a smooth underlying spacetime. The wave 
amplitude is 



M + (M - a z ) 



2\l/2 



(2) 



GAL 



quad 



re- 



(3) 



The area of the event horizon is 8irMr+. 



where Q is the quadrupole moment of the source, r is 
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FIG. 1 Lines of force for plane gravitational waves propagat- 
ing along the z axis. The wave on the left is purely in the 
+ polarization state, and the one on the right is purely in 
the x polarization state. The gravitational waves produce 
tidal f orces in planes transverse to the propagation direction. 
From (HE ramovici et all \l99$) . Reprinted with permission 
from A A AS. 



the distance from the source, and M qua( j is the mass in 
the source that is undergoing quadrupolar changes. This 
shows that the strongest gravitational waves will be pro- 
duced by large masses moving at high velocities, such as 
binaries of compact stars and black holes. 

A gravitational wave is purely transverse, acting tidally 
in directions perpendicular to its propagation direction. 
When a gravitational wave impinges on a detector of 
length scale L, it produces a length change in that de- 
tector SL/L ~ h/2. By substituting in typical values for 
compact objects in the Universe into Eq. ©, one can 
see that astrophysical sources typically yield wave ampli- 
tudes of h < 10 -21 at the Earth. Consequently, precision 
measurements are needed to make detections. 

Gravitational waves have two polarization compo- 
nents, known as h + and h x for a linearly polarized 
wave. Figure Q] shows the corresponding lines of force 
for sinusoidal gravitational waves propagating along the 
z axis. In the purely + polarization on the left, the wave 
stretches along one axis and squeezes along the other, 
alternating sinusoidally as the wave passes. The x po- 
larization wave on the right acts similarly, stretching and 
squeezing along axes rotated by 45°. In general, a gravi- 
tational wave is a superposition of these two states, con- 
veniently written as a complex waveform strain h, where 



(4) 



B. Astrophysical Black Holes 



Astronomers have found evidence for black holes 
throughout the Universe on a remarkable range of scales. 
The smallest of these, stellar black holes, have masses in 
the range ~ (3 — 3O)il/ and form as the end-products 
of massive star evolution. There is good observational 
evidence for the existence of stellar black holes, based on 
dynamical measurements of the masses of compact ob- 
jects in transient systems that undergo X-ray outbursts. 
Since neutron stars cannot have masses > 3M Q , com- 
pact objects more massive than th is must be black holes 
( Remillard and McClintockLl2006h . 



Intermediate-mass black holes (IMBHs) have masses 
in the range - (10 2 - 10 4 )M o . IMBHs may 
form as the result of multiple mergers of smaller 
objects in the centers of dense stellar cl usters 
in the present Universe (iMiller and Colbertl . 120041 : 
iPortegies Zwart and McMillan! . 120021) . assuming mass 
loss f rom stellar winds is not significant ( Glebbeek et all 
l2009h . They may also arise from the evolution of very 
massive stars early in the history of the Universe, form- 
ing black-hole "seeds" in the centers of massive ha- 
los (the precursors of the galaxies we see today) early 
i n the history of the U niverse, to redshifts z > 10 
( Madau and Reed . l200ll ). Currently the best observa- 
tional evidence for IMB Hs comes from models of ultra- 
luminous X-ray sources (IColbert and Milled . 120051) . 

Finally, massive black holes (MBHs) have masses in 
the range ~ (10 4 — 10 9 )A/ Q and are found at the cen- 
ters of galaxies, including our own Milky Way galaxy. 
The observational case for the existence of MBHs is quite 
strong, based on dynamical models of stars and gas be- 
lieved to be moving in the potential well of the central 
MBH dDesroches et all. l2009t iFerrarese and Ford|. 120051: 



iKormendv and Richstonel ." 1995t iRichstone et all 19981 ). 

Black-hole binaries are binary systems in which each 
component is a black hole. As mentioned above, we focus 
here on comparable-mass binaries, which are expected to 
produce the strongest gravitational-wave signals. Stel- 
lar black-hole binaries may form a s the result of binaries 
compo sed of two massive stars; see lBulik and Belczvnskil 
(120091) and references therein. Stellar black-hole bina- 
ries may also arise from dynamical processes in which a 
black hole is captured into an orbit around another black 
hole i n dense stellar e nviro nments ( Miller and Lauburd . 
120091 ; lO'Learv et all I2007D . IMBH binaries can also 
form through dyna mical processes in stellar clusters 
( Gii rkan et al ., 200G) and from mergers of massive halos 
at high redshifts. Since both stellar and IMBH binaries 
are "dark" - that is, they are generally not surrounded 
by gas which might produce electromagnetic radiation - 
we have few observational constraints on these types of 
black-hole binaries. This situation will change dramat- 
ically, however, with the d etection of gravitational ra- 
diation from these systems ( Bulik and Belczvnskil . l2009t 
iMillerl 120091) . as gravitational waves bring direct informa- 
tion about the dynamical behavior of the orbiting masses 
and do not rely on electromagnetic emissions from nearby 
matter. 

Since essentially all galaxies are believed to contain an 
MBH at the center and to undergo a merger with an- 
other galaxy at least once during the history of the Uni- 
verse, MBH binaries can aris e when t heir host galaxies 
merge ( Begelman et all\l98(t ); see also iDiorgovski et al\ 
(2008) and references therein. However, due to the vast 
cosmic distances involved, and the small angular sepa- 
rations on the sky expected for MBH binaries, only a 
few candidates are currently known through electromag- 
netic observations (iKomossa . 120031 ; iKomossa et al. 1. 120031: 
iRodriguez et all 120061 ). When they form, MBH bina- 
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ries typically have relatively wide separations, and the 
gravitational radiation they emit is very weak and insuf- 
ficient to cause the binary to coalesce within the age of 
the Universe. However, various processes such as gaseous 
dissipation and iV-body interactions with stars can re- 
move orbital energy fro m the binary and cause the black 
holes to spiral togethe r ( Armitage and Nataraianl 120021 : 
IGould and RH l2000l): see also iBerentzen et olT(|2009| i 
and Colpi et all (|2009T) and references therein. Eventu- 



19571 : iTeukolskvl Il973t IZerillil . 1 1970ft form the basis 
of ringdown calculations, producing gravitational wave- 
forms in the shap e of exponentially damped sinusoids 
(jBerti et a/.l . l2009t iLeaveH Il986ft . 

The characteristic gravitational-wave frequency of a 
quasicircular black-hole binary, produced by the domi- 
nant (highest-order) quadrupolc component, is 



/gw ~ 2/ orb ~ (M/R 



3\l/2 



(5) 



ally, the black holes reach separations at which gravita- 
tional radiation reaction becomes the dominant energy- 
loss mechanism, leading to the final coalescence of the 
black holes and the em ission of strong gravitational waves 
(jSesana et all l2009bft . 



C. Gravitational Waves from Black-Hole Binaries 

Mergers of comparable- mass black- hole binaries are ex- 
pected to be among the strongest sources of gravitational 
waves. This final death spiral of a black- hole binary en- 
compasses three stages : insp i ral, merger, a nd ringdown 
(|Flanaean and Hughes! . Il998t lHughesl . l2009ft . 

In the early stages of the inspiral, the orbits of most 
astrophysical black- hole binaries will circ ularize due to 
the emission of gravitatio nal radiation ( Petersl . Il964t 
iPeters and Mathewsl . Il963ft . During the inspiral, the or- 
bital time scale is much shorter than the time scale on 
which the orbital parameters change; consequently, the 
black holes spiral together on quasicircular orbits. Since 
the black holes have wide separations, they can be treated 
as point particles. The inspiral dynamics and waveforms 
can be calculated using post-Newtonian (PN) equations, 
which result from a systematic expansion of the full Ein- 
stein equations in powers of \e^_jPj^j^_GM/Rc 2 , where m HISTORICAL OVERVIEW 



where / or b is the orbital frequency. Astrophysical black- 
hole binaries produce gravitational waves that span three 
frequency regimes, dependin g on the black-hole masses 
( Flanagan and Hughesl |2005[) . Stellar black- hole binaries 
and the lower mass end of the IMBH binaries radiate in 
the high frequency band, /gw ~ (10 — 10 4 )Hz, which is 
already being observed by gro und-based laser interf erom- 
eter detectors such as LIGO ( Abbott et all [2009b[ ) . and 
will be obser ved by the advanced detectors by ~ 2016 
( Smith! . 120091 ) . Low frequency gravitational waves cover 
the band /gw ~ (10~ 5 — l)Hz and will be observed by 
the space-based la ser interferomet er LISA, currently un- 
der development (jjennrichl . |2009() . MBH binaries with 
masses M ~ (10 4,5 — 10 7 )Af Q will be very strong sources 
for LISA, wit h the lower mass sy stems visible out to red- 
shifts z > 10 (jArun et all l2009al ) ; the inspiral s of IMBH 
binaries will also be detectable ( Millerl . 12009ft . Finally, 
the very low frequency band /gw ~ (1 0~ 9 — 10~ 7 )Hz 
will b e observed by pulsar timing arrays (jVerbiest et all 
12009ft . This band is expected to dominated by gravita- 
tional waves fr om a very large pop ulation of unresolved 



MBH binaries (Sesana et al 



discrete sources (iScsana et 



20081) w ith possibly a few 



2009aft . 



R is the binary separation ( Blanchetl 12006ft . The inspiral 
phase produces gravitational waves in the characteristic 
form of a chirp, which is a sinusoid with both frequency 
and amplitude increasing with time. 

As the black holes spiral inward, they eventually reach 
the strong-field, dynamical regime of general relativity. 
In this merger stage, the orbital evolution is no longer 
quasi-adiabatic; rather, the black holes plunge together 
and coalesce into a single, highly distorted remnant black 
hole, surrounded by a common horizon. Since the point- 
particle and PN approximations break down, numerical 
relativity simulations of the Einstein equations in three 
dimensions are needed to calculate the merger. Due to 
the difficulty of these simulations, the resulting gravi- 
tational waveforms were completely unknown until re- 
cently. 

Finally, the highly distorted remnant black hole settles 
down into a quiescent rotating Kerr black hole by shed- 
ding its nonaxisymmetric modes through gravitational- 
wave emission. We call this process the "ringdown," 
in analogy to how a bell that has been struck sheds its 
distortions as sound waves. Variou s analytic techniques 
of black- hole perturbation theory ( Regge and Wheelerl . 



The quest to calculate the gravitational-wave signals 
from the merger of two black holes spans more than four 
decades and encompasses key developments in theoreti- 
cal and experimental general relativity, astrophysics, and 
computational science. In this section, we begin by de- 
lineating these threads in general terms, and then turn 
to a more detailed account of select milestones along the 
path toward successful simulations of black-hole mergers. 



A. Setting the Stage 



At th e end of the 18th century iMichelll (| 17841 ) and 
lLaplacel (1796) first speculated, using Newtonian grav- 
ity, that a star could become so compact that the es- 
cape velocity from its surface would exceed the speed of 
light. In the 20th century, scientists realized that such 
black holes could form as the final state of total gravita- 
tional collapse in general r e lativit y ( Harrison et all [19651: 
lOppenheimer and Snvderi . Il939ft : John Wheeler would 
later po pularize the term "b l ack hole" to describe s uch a n 
object (jMisner et all 120091: iRuffini and Wheeled . Il97ll) . 
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Beginning in the 1960s, many highly energetic astrophys- 
ical phenomena were discovered with physical properties 
pointing to extremely strong gravitational fields as un- 
derlying mechanisms; among these are quasars and X-ray 
binaries such as Cygnus X-l ( Overbeck and TananbaurrJ . 
119681 : lOverbeck et all Il967h . the first credible black-hole 
candidate. As discussed in Sec. III.B1 today astrophysical 
black holes arc believed to exist on a vast range of scales 
throughout the Universe, and black-hole binaries are con- 
sidered to be strong sources of gravitational waves. 

Einstein's equations of general relativity form a cou- 
pled set of nonlinear partial differential equations, in 
which dynamic curved spacctime takes the role of New- 
ton's gravitational field and interacts nonlinearly with 
massive bodies. Gravitational waves were first recognized 
as solutions to the linearized, weak-field Einstein equa- 
tions early in the past century. By mid-century, gravi- 
tational waves were recognized as real physical phenom- 
ena, carrying energy and being capable of producing a 
response when impinging on a detector. This develop- 
ment spawned a major branch of experimental general 
relativity, with concepts for the first gravitational- wave 
detect ors appearing in the 1960s; see lCamp and Cornish! 
(|2004l ) for a review. The discovery of two neu tron stars 
in a binary system by iHulse and Tavlorl ( 19751 ) provided 
an astrophysical laboratory for the first indirect detec- 
tion of gravitational radiation. Decades of observa- 
tion revealed the binary orbit to be shrinking by pre- 
cisely the amount expected if the system were emit- 
ting gravh^ationaI_waves accor ding to general relativ- 
ity ( Weisberg and Tavlorl 120051 ); Hulse and Taylor were 
awarded the Nobel Prize in 1993. Today, the progno- 
sis for direct detection of gravitational waves is excellent, 
with the first events expected from the advanced ground- 
based interferometers around the middle of this coming 
decade. 

As you might expect, Einstein's equations pose 
formidable obstacles to anyone who would dare to probe 
the physics within. Throughout most of the 20th cen- 
tury, relativists uncovered a fairly small number of exact 
solutions by exploiting symmetries, and made progress 
toward more general problems using various perturba- 
tive expansions. By the 1960s, computers had become 
powerful enough to encourage attempts at solving Ein- 
stein's equations numerically, to uncover physics beyond 
the realm of perturbation theory. The subsequent devel- 
opment of numerical relativity has been made possible 
in part by continued increases in computer power and 
advances in algorithms and computational methods. 

Most numerical-relativity simulations start with the 
idea of decomposing four-dimensional spacctime into 
a stack of curved three-dimensional spacel ike slices 
threa ded by a congruence of timelike curves ( York Jrl 
Il979l ). Arnowit, Deser, and Misner (ADM) pioneered this 
"3+1" approach as the basis for a cano nical formulation 
of th e dynamics of general relativity (jArnowitt et all . 
11962ft . In this Hamiltonian formulation, the three-metric 
-fij on the spatial slices takes the role of the "configura- 



tion variables." Quantities based on the extrinsic curva- 
ture Kij , which is roughly the time derivative of 7^ , play 
the role of "conjugate momenta." Variation of an action 
with respect to 7^ produces a set of six first-order evo- 
lution equations for the conjugate momenta; varying the 
momenta gives six first-order evolution equations for 7^' . 
ADM also introduced four Lagrange multipliers as freely 
specifiable gauge or coordinate conditions, representing 
the four coordinate degrees of freedom in general relativ- 
ity. Variation of these Lagrange multipliers yields four 
equations that must hold on each slice: a Hamiltonian 
constraint, and three momentum constraints. 

Originally intended as a tool for quantizing gravity, the 
ADM formalism later became the basis for most work in 
numerical relativity. As we discuss in Sec. IIVI below, 
key elements in this approach are solving the Cauchy 
problem, beginning with the initial data on a spacelike 
slice, and then evolving that data forward in time. The 
constraint equations form the basis for this initial value 
problem. Appropriate choices for the gauge conditions 
are crucial ingredients for today's successful black-hole 
merger simulations. 



B. Numerical Relativity Milestones 



lHahn and Lindquistl (|l964l made the first known at- 
tempt to simulate the head-on collision of two equal- 
mass black holes on a computer in 1964, using a two- 
dimensional axisymmetric approach. Their simulation 
ran for 50 timesteps to a duration of ~ 1.8M; at 
this point, they decided that the simulation was no 
longer accurate enough to warrant continued e volution , 
and s topped the code. Smarr and Eppley ( Epplevl . 
119751 : ISmarrl . Il975l . Il977l : ISmarr et all . 119761) returned 
to this problem in the mid-1970s, again employing two- 
dimensional axisymmetry but now using the ADM for- 
malism, specialized coordinates, and improved coordi- 
nate conditions. Although they encountered problems 
with instabilities and large numerical errors, they man- 
aged to evolve the collision and extract information about 
the emitted gravitational waves. Smarr and Eppley had 
used the most powerful computers of their day. Going to 
the next step, orbiting black holes in three dimensions, 
was deemed to be not feasible at the time, due to un- 
resolved questions about the instabilities and insufficient 
computer power. Consequently, the black- hole merger 
problem lay dormant for over a decade. 

In the 1990s, attention returned to black-hole mergers 
as the LIGO project began to move forward, and black- 
hole mergers were recognized as the strongest sources 
for this detector. Since the signal-to-noise ratio for 
such ground-based detectors is fairly modest, having a 
template for the merger waveform is a key part of the 
data analysis strategy. Numerical relativists revisited the 
problem of colliding two black holes head- on with modern 
techn i ques and more p o werfu l computers ( Anninos et all , 
Il993l : iBernstein et all . 1 1994) . In the mid-1990s, the 
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National Science Foundation funded a Computational 
Grand Challenge grant for a large multi-institution col- 
laboration aimed at evolving black-hole mergers in three 
dimensions and calculating the resulting gravitational 
waveforms. Around the same time, a large and very ac- 
tive numerical relativity group arose at the newly-formed 
Albert-Einstein Institut (AEI) in Potsdam, Germany. 
During the late 1990s and into the early 2000s, two de- 
velopments on the experimental side further increased 
the desire for black-hole merger simulations: the ground- 
based gravitational-wave detectors started taking data, 
and interest grew in LISA and its potential for observ- 
ing gravitational waves from massive black-hole binary 
mergers. 

While no one expected the task at hand - develop- 
ing computer codes to solve the full Einstein equations 
in three dimensions for the final few orbits and merger 
of two black holes - to be simple, numerical relativists 
found that the problem was far more difficult than antic- 
ipated. Producing a waveform useful for gravitational- 
wave detection purposes typically would require running 
a simulation for a duration of several hundred M . How- 
ever, a variety of instabilities plagued the codes, causing 
them to crash well before any significant portion of an 
orbit could be achieved. 

Nevertheless, during a period encompassing a little 
over a decade, much important work was accomplished 
that laid the foundations for later success. Key mile- 
stones include these developments: 

• initial data for b in ary b l ack h ol es ne a r the 
ISCO (iBaumgartd. 120001 ICookL 11994 120021: 
ICook and Pfeifferl |2004) ; 



• new methods for representing black holes 
on computational gr ids suc h as punctures 

excision 



( Brandt and Briigmannl. Il997ft and 
( Alcubierre and Briigmann, 1200 it Anninos et al. . 




Seidel and Suenl . 1992t Shoemaker et al . 



• recognition of the importance of hyperbol 
icity in formulating the Einstein equations 



lcity in iormuiatmg t nc tLmstem eg r 
for numerical s olution (lAbrahams et all , 
Anderson et all 1997 ; Bona and Massol 



Friedrich and Rendalll . boOOl ) 



1997 



1992: 



• improved formulations of the Einstein equations 
i 1 l I 1 I 1 I 

(Baumgartc and Shapiro, 1998 : | Nakamura et all 
119871: IShibata and Nakamural Il995l ); 



• fully three-dimensional evolution c odes and their 
use in evolving distorted black hol es (jBrandt et al\ . 
120031: ICamarda and Seidell Il999t) . boosted black 
holes (ICook et all | |998l ). head-on collisions 
oobTT a nd grazing c ol lisions 
l2001at iBrandt et all . l200Ct 



( Soer hake et all 
(I Alcubierre et all 
lBrikmannl ll999) 



• coordinate conditions that keep the slices from 
crashing into singularities and the spatial coordi- 
nates from fall ing into the black holes as the evolu- 
tion proceeds ( Alcubierre et all l2003t iBona et all , 
[1991); 

• the Cactus Computational Toolkit 1 , which pro- 
vided a framework for developing numerical- 
relativity codes and analysis tools used by many 
groups; 



moaern aaapuve mesn renncmenr nmic-amerenc 
(inclu ding Carpet 2 wi th Cactus, BAM rtBrugmanr 
19991) . and Paramesh dMacNeice et al.). 120001) 1. am 



• modern adaptive mesh refinement fini te-difference 

and 

multi-domain spectral (jPfeiffer et all l2003f ) infras- 
tructures for numerical relativity. 

Throughout this period, the length of black-hole evo- 
lutions gradually increased. Improvements in the for- 
malisms allowed simulations of single black holes and, 
later, two black holes to increase in duration to > 10M. 
The addition of new slicing and shift conditions again 
increased the evolution times to > 30M. The Lazarus 
project took a novel approach to combine these rela- 
tively short-duration binary simulations with perturba- 
tion techniques for the late-time behavio r to produce 
a hy b rid model fo r a black-hole mer ger dBaker et all 
I2000L 120011 l2002allbh . By late 2003. iBrugmann et al\ 
(|2004l ) carried out the first complete orbit of two equal- 
mass, nonspinning black holes. While this simulation 
lasted ~ 100M, the code crashed shortly after the or- 
bit was completed and the gravitational waves were not 
extracted. Overall, progress was slow, difficult, and in- 
cremental. However, the situation was about to change 
dramatically. 



C. Breakthroughs and the Gold Rush 

In early 2005, Frans Pretorius electrified the relativ- 
ity community when he achieved the first evolution of 
an equal-mass black-hole binary through its final or- 
bit, merger, and ringdown using techniques very diffcr- 
ent from those em ployed by the rest of the community 
( Pretoriusl l2005al ) . Later in 2005 the groups at the Uni- 
versity of Texas at Brownsville (UTB) and NASA's God- 
dard Space Flight Center independently and simultane- 
ously discovered a new method, called "moving punc- 
tur es," that also produced successful black-hol e merg- 
ers (|Baker et aUl2006btlCampanelli et aZl . l2006al ). Their 
presentations were given back-to-back at a workshop on 
numerical relativity, to the amazement of each other and 
the assembled participants. 

Since the moving-puncture approach was based on un- 
derlying techniques used by several other groups, it was 



http: //www. cactuscode . org/ 



http: //www. carpetcode . org/ 
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rapidly and readily adopted by most in the community, 
producing a growing avalanche of results. 2006 was a year 
of many firsts: the first simulations of unequal-mass black 
holes and th e accom panying recoil of the remnant hole 
( Baker et all 2006c ), the first m ergers of spinning black 
holes ([Campanelli et all 2006df). t he first long waveforms 



(~ 7 orbits) ([Baker et all l2007d) . the first comparisons 



with PN results ([Baker et ai , 2007dt iBuonanno et all 
l2007al ). and the fir st systematic pa r ameter study in nu- 
merical relativity ( Gonzalez et all l2007b|) . The year 
2007 opened with the discovery of "superkicks" - re- 
coil velocities exceeding 1000 km s _1 ( Campanelli et ali , 
l2007bl: iGonzalez et aLl . l2007ah . and in 2008 a black-hole 
binary merger w i th ma ss ratio q = 10 was accom p lished 
([Gonzalez et all |2009( ). while ICampanelli et all ([2009h 



carried out the first long-term evolution of generic spin- 
ning and precessing black-hole binaries. As this article 
was being written in 2009 the state of the art contin- 
ues to advance, with the first simulations o f black-hole 
mergers using spec tral numerical techniq ues (|Chu et all 
l2009t IScheel et all 120091: ISzilagyi et all , |2009j), and the 
first steps towards m odeling the flows of gas around the 
merging black holes ( van Meter et all l2010bl ) . Applica- 
tions of the merger results in areas such as comparisons 
with PN expressions for the waveforms, astrophysical 
computations of black- hole merger rates, and the devel- 
opment of templates for gravitational- wave data analysis 
have accompanied these technical developments in the 
simulations. The study of black-hole mergers using nu- 
merical relativity is thriving indeed. 



IV. NUMERICAL DEVELOPMENT 

In this section, wc discuss the mathematical and nu- 
merical foundations underlying current black-hole merger 
simulations, highlighting the key issues involved in 
achieving successful evolutions. For more det ailed treat- 
ments , w e direct the interes ted reader to lAlcubierre] 
(120081) or (lGourgoulhonll2007D . 



The metric characterizes the geometry of spacetime by 
giving the infinitesimal spacetime interval ds through the 
following definition: 



ds 2 = g^dx^dx" , 



(7) 



where we use the Einstein summation convention which 
implies summation over all values of a given index that 
appears twice in an expression. Many physical impli- 
cations of the metric are immediately apparent from 
Eq. ([7|). For example, when ds 2 = 0, the resulting metric- 
determined relationship between the time and space co- 
ordinates yields the paths that light rays must follow in 
this spacetime. 

The dependence of Einstein's tensor on the metric 
can be simply illustrated in coordinates known as "har- 
monic" ; in vacuum Eq. ^ takes the form: 



a 9fi,u - t^u = 0, 



(8) 



where □ is the flat-space wave-operator and i M „ repre- 
sents all terms nonlinear in the metric. If t^ is inter- 
preted as an effective source term, this is a simple wave 
equation. The familiar form of this equation suggests 
that its Cauchy problem can be solved by specifying the 
metric and its first derivative on an initial, spatial sur- 
face and then integrating in time, as for an ordinary wave 
equation. 

Einstein's equations admit various formulations and 
coordinate conditions, which should be tailored to the 
problem at hand - in this case, numerical simulation 
of black-hole spacetimes. Regardless of these choices, 
current numerical practices universally involve an initial 
three-dimensional slice of spacetime that is evolved for- 
ward in time. Here, we review the history and current 
most common choices of initial data, black-hole repre- 
sentations, formulation, coordinate conditions, and some 
details of the numerics. The impetus, and most impor- 
tant outcome, of all these developments is the ability to 
generate gravitational waveforms from black-hole binary 
sources over many cycles. 



A. Einstein's Equations 

The central task of numerical relativity is solving Ein- 
stein's field equations 

G>j, = SttT^, (6) 

where Einstein's tensor G^ represents the curvature of 
spacetime, the energy-momentum tensor T^ v contains 
the matter sources, and fj,, v = 0,1,2,3. By conven- 
tion, an index /x = selects a " time" component, and 
fj, = 1, 2, 3 selects a "space" component. G^ depends on 
the first and second derivatives of the metric tensor g^. 
For vacuum black-hole spacetimes, T^ v = 0. Note that 
all of the tensor fields discussed here are symmetric in 
the indices, e.g. g^ v = g vlll . 



B. The Cauchy Problem 

The Cauchy problem concerns solution of the field 
equations given initial data specified on an initial (typi- 
cally spatial) hypersurface. In practice, the Cauchy prob- 
lem is more conveniently investigated in a "3+1" formu- 
lation, explicitly based on a foliation of the spacetime 
into three-dimensional spatial slices parametrized by a 
time coor dinate. A common 3 +1 formulation inspired 
by ADM (lArnowitt et all Il962f ) divides up the compo- 
nents of the metric according to their relationships with 
space and time, such that the line element takes the form: 

ds 2 = {-a 2 + ftp^dt 2 + 2l3 i dtdx i + - ,//.-•'(/.,•'. (9) 

where a is called the lapse function, (3 l the shift vector, 
and 7ij = is the spatial three-metric. We write the 
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time coordinate x° — t and the spatial coordinate indices 
i, j = 1, 2, 3. Note that contraction with g^ v or its inverse 
g^ v is used to lower or raise indices of four-dimensional 
tensors, respectively, while 7^ or its inverse r ) lJ is used to 
lower or raise indices of three-dimensional tensors. The 
lapse and shift represent coordinate freedom in the met- 
ric; we can choose these quantities arbitrarily. However, 
since the three-metric 7^ (and its first and second spa- 
tial derivatives) determines the intrinsic curvature of the 
slice, it carries the information about the gravitational 
field and thus is constrained by the physics. 

The meaning of the lapse and shift can be understood 
by considering two successive spatial slices separated by 
an infinitesimal time interval dt (Fig. [5]). An observer 
along a vector normal to the first slice will measure an 
elapsed proper time of dr — adt in evolving to the second 
slice, and a change in spatial coordinate of dx % — ~(3 l dt. 

Since Einstein's equations are second order in time, we 
must also specify the initial time derivative of the three- 
metric. Rather than specifying this derivative directly, 
we define a new quantity: 



K i:j = (dtlij 



(10) 



where dt = d/dt is an ordinary partial derivative and Di 
is the spatial covariant derivative. Note that the space- 
time covariant derivative V M is a partial derivative with a 
correction such that it transforms as a vector and satisfies 
VaS^i/ = 0. Di is the projection of onto the spatial 
slice and is equivalent to a three-dimensional covariant 
derivative formed from 7^ . For the case of Euclidean 
normal coordinates, a = 1 and f3 l = 0, Eq. (fTU|) reduces 
to the simple expression Kij = — (1/2)9* 7^. 

If we define a unit vector normal to the spatial slice, 
we can show that Eq. (|10[) is equivalent to = —DiUj. 
As suggested by this expression, is a measure of the 
change of the normal vector as it is transported along the 
slice. In this way gives an extrinsic measure of the 
curvature of a three-dimensional spatial slice with respect 
to its embedding in four-dimensional spacetime. It is 
therefore called the extrinsic curvature. Depending on 
the formulation, the extrinsic curvature might or might 
not come into the evolution equations; however, is 
almost universally utilized when calculating initial data. 

Of the ten component equations of Eq. © , six deter- 
mine the time-evolution of the metric, while four must 
be satisfied on a spatial slice at any given time, and are 
thus constraint equations. With an appropriate choice of 
time coordinate, and assuming vacuum spacetime, these 
four constraint equations are equivalent to the condition 
Gov = 0. The time-time component Goo = is called the 
Hamiltonian constraint, and the time-space components 
Goi = are called the momentum constraint. These take 
the form of conditions on the extrinsic curvature: 



(3) 



adt 



^jif-J^dt^ 


Jx" 

\ # — 


1 

r ' 







t = dt 



t = 



R + K - 
Dj (K i3 



KijK 13 = 



(11) 
(12) 



FIG. 2 The 3+1 split into space and time. Two spatial slices 
at t = and t = dt are depicted, a is the lapse, and adt 
represents the proper time lapse between slices. /3 a is the 
shift, and fi a dt represents the amount by which the spatial 
coordinates shift between slices. n a is normal to the slice at 
t = 0. If a ray parallel to n a intersects the t = slice at a 
point x a , then it will intersect the t = dt slice at x a — /3 a dt. 
t a = an a + /3 a is a coordinate time vector. If a ray parallel 
to t a intersects the t = slice at point x a , then it will also 
intersect the t — dt slice at x a . 



where is the three-dimensional Ricci scalar associ- 
ated with the three- metric 7^, and K = traceifjj = 
7 y , sometimes called the mean curvature. These con- 
straint equations must be solved in order to obtain an 
initial spatial slice consistent with Einstein's equations. 



C. Representing Black Holes in Numerical Spacetimes 

How does one represent an exotic object such as a 
black hole in a numerical simulation? In particular, how 
can one use finite computational methods repeatedly to 
model an object which, analytically, contains physical 
and/or coordinate singularities? Fortunately, two suc- 
cessful strategies have emerged to meet this challenge. 

The unusual topology of bla ck holes offers one way out. 
As lEinstein and Rosen ( 1935| ) originally showed, a black 
hole can be considered a "bridge" or "wormhole" con- 
necting one Universe or "worldsheet" , to a second world- 
sheet (see Fig. [3]). Exploiting this topology, a continuous 
spatial slice which avoids the physical singularity con- 
tained within the event horizon of each black hole can 
be constructed as follows. Starting with a spatial slice 
of Schwarzschild spacetime, remove the interior of the 
event horizon. Identify the resulting spherical bound- 
ary with the spherical boundary of an identical copy of 
this space. The two-dimensional analog would be to take 
two sheets of paper, cut out a disk from each, and then 
glue the resulting circular edges together. Each sheet, 
or copy of Schwarzschild spacetime in this example, is 
called a worldsheet. The identified spherical boundaries 
connecting the worldsheets form what is referred to as 
the "throat" of the wormhole. 

To complete this construction, we require an appropri- 
ate coordinate system to continuously cover both world- 
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sheets. Brill and Lindquist discovered coordinates that 
will prove convenient, in which the three-metric of a par- 
tic ular spatial slice of Schwar zschild spacetime is given 
by (jBrill and Lindouistl . fl963l) 



7tf = 1 



(13) 



where r is a radial coordinate. In these coordinates, the 
event horizon is at r = m/2. We can consider each of 
the worldsheets described above as being separately la- 
beled by such coordinates. Designate one worldsheet as 
A and the other as B, and call their radial coordinates r a 
and Tb respectively. Since the interior of the event hori- 
zon has been removed from each worldsheet, assume that 
Ta > m/2 and tb > m/2. The metric on each worldsheet 
has the form of Eq. (fTS]) . As noted by Brill and Lindquist, 
this form of the metric is unchanged by the transforma- 

when r = m/2. So if we 



tion r' = m 2 /Ar, and r' 
define a new coordinate r' by 



ta 



r tor r > —, r B 



2? tmr - r 



(14) 



mapping spatial infinity on worldsheet B to r' — 0, we 
obtain a single continuous coordinate system that applies 
to worldsheet A for r' > m/2 and worldsheet B for r' < 
m/2. 

This avoidance of the physical singularity comes at 
the expense of a coordinate singularity in the metric 
at r' = 0, called a "puncture". However, it turns out 
this coordinate singularity can be confined to a single 
scalar variable, as suggested in Eq. (p~3|) . With a suitable 
change of variables and other means, numerical simu- 
lations haven proven capable of handling this irregular 
scalar field. Thus the "puncture" method is one way to 
represent a black hole that is amenable to computation. 

Another strateg y is excision, first proposed by Unruh 
( Thornburgl . Il993h . Given that no physical information 
can escape an event horizon to influence the exterior, 
the interior of a black hole can in principle be excised 
from the computational domain 3 . This relies on the fact 
that all physical information propagates inwards from the 
event horizon towards the physical singularity, i.e. light- 
cones tilt inwards, and nothing physical propagates out- 
ward from the horizon. Extrapolation can be used for the 
boundary condition of the excised region, and any non- 
physical numerical error that escapes the horizon should 
be negligible. This approach is not constrained to the 
particular coordinates required by the puncture method, 
but it can be more difficult due to the need for precise 
positioning of the excision boundary and accurate ex- 
trapolation. Excision was used success fully for orbiting 
binary simulations bv lPretor ius (20053) and continues to 



3 As unphysical coordinate "gauge modes" may couple to physical 
modes, we also assume no superluminal coordinate effects are 
present. 




- throat r' = r/, 



FIG. 3 In the wormhole representation of a black hole, the 
initial slice typically just touches the horizon. The upper 
"sheet" represents the exterior space, while the lower sheet is 
a duplicate, joined to the upper sheet by a "throat". 



be us ed by the Ca l tech-Cornel l dual-coordinate spectral 
code (IScheel et ali . 120091 . 12006T) . 

Either of these methods can be useful in representing 
black holes on the initial data slice. Surprisingly, we will 
also see that either of these representations can be made 
robust enough to persist as the black holes evolve. 



D. Initial Data 

The starting point for successful simulations of black- 
hole mergers is finding initial data for astrophysically re- 
alistic inspiralling black holes. If we were simulating the 
orbits of stars in Newtonian gravity, this would be a sim- 
ple procedure. For example, we could simply specify the 
masses and spins of the stars, along with their positions 
and velocities on orbits derived from the dynamics of 
point particles, and then evolve the system numerically 
to allow it to "relax" into orbits appropriate for bodies of 
finite size. In general relativity, however, the initial data 
must satisfy the constraint equations ([TTj) - (fl"2"|) . which 
are cast in terms of the 3-metric 7^ and the extrinsic 
curvature Kij. Since there is no obvious, or natural, 
connection between these field variables and the astro- 
physical properties of inspiralling black holes, obtaining 
suitable initial data is a major challenge. 

Building on the earlier work of iLichnerowiczl ( 19441 ) . 
York developed a general procedure for solving the con- 
straint equations to prod uce initial data for the Cauchy 
problem in the 1970s (see I York Jr.1 ( 19791 ) for a review). 
This approach generally requires solving a coupled el- 
liptic system of four nonlinear field equations. We can 
break this problem down into more manageable pieces 
with some simplifying assumptions. While these choices 
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do come at the cost of some loss of generality and astro- 
physical realism in the initial data for two black holes, 
it has been seen that for sufficiently long evolutions the 
final orbits and waveform signatures from the black-hole 
evolution are largely insensitive to this level of detail in 
the initial data, at least in the case of equal-mass non- 
spinning black holes. 

The first simplification is to choose traceless extrinsic 
curvature, K = 0. With this, the Hamiltonian constraint 
is decoupled from the momentum constraints, and can 
be solved separately. To find solutions corresponding to 
multiple black holes, we generally further assume that 
the initial slice is conformally flat. That is, the three- 
metric is the product of a scalar conformal factor with a 
flat metric, 



1 >j 



1>% 



(15) 



With this, the problem reduces to solving a (typically 
n onlinear) equation fo r the s calar field ip. 

iBrill and Lindouistl (jl963h found a simple solution, 
representing N black holes momentarily at rest, which 
gives 



N 

i+y^i 



(16) 



where is the mass associated with the ith black hole 
and Ti represents its coordinate center. Each r, corre- 
sponds to the location of a puncture, as described in the 
last section. This is a valuable solution, but not gen- 
erally useful, because these black holes lack momentum 
and spin. 

An explicit solution by Bowen for the momentum con- 
straint was the last crucial step in defining a procedure for 
calculating initial data for multiple black holes with spec- 
ified linear and angular momenta ( Bowen and York Jr. I . 
1 1 980h . The Bowen- York prescription employed a two- 
sheet ed topol o gy fou nd by Misner for the black-hole inte- 
riors ()Misnerl . fl963h . Later Brandt and Brugmann gen- 



eralized the procedure f or the Brill-Lindquist topology 
( Brandt and Brugmannl . Il997l ). The Brandt-Briigmann 
puncture data for arbitrary momenta and spin is now 
widely used because of its ease of implementation. 

More recently York developed another modeling 
ansatz, kn own as the Conforma l -Thin - Sandwich (CTS ) 
approach dPfeiffer and York Jrl 120051: lYork Jrl Il999h . 
which has certain additional advantages. For example, 
it turns out the conventional Brandt-Briigmann punc- 
ture ca nnot yield a spin parameter (a/M) greater than 
0.93 (iDain et all 120021) . while CTS data can go higher 



( Lovelace et all 2008[ ). As discussed above, the spatial 



metric is chosen to be conformally flat, then instead of 
providing an ansatz for the extrinsic curvature Kij, the 
initial time-derivative of the conformal metric is speci- 
fied (generally to vanish). In addition, a condition can 
be imposed on the slicing (see below). The result is a 
coupled system of elliptic equations which is solved to 



enforce the constraints and optionally the choice of slic- 
ing condition. Boundary conditions may additionally be 
supplied to enforce rotational motion. 

Either approach provides an ansatz for constructing 
three-dimensional binary black-hole initial field data for 
a specified choice of particle-like parameters, masses, po- 
sitions, momenta, and spins. Inevitably there are differ- 
ences from the nearly quiescent evolving systems that we 
seek to represent. Generally there is some level of spu- 
rious radiation generated from a period of initial tran- 
sient dynamics through which the system relaxes to be- 
come quiescent on sub-orbital time scales. In particu- 
lar, extraneous radiation content is an unavoidable conse- 
quence of conformal flatness, which post-Newtonian anal- 
ysis has shown must deviate from the physi c ally r elevant 
inspiralling binary solution ( Damour et all |2000( ) . This 
spurious radiation is seen in plots of the gravitational 
waveforms produced by mergers; see Fig. [9] and other 
waveform plots in Sec. IV.B.2I Often the simulations will 
also undergo a period of initial gauge-evolution which, 
though physically inconsequential, may affect the quality 
of the simulation num erically. For puncture initial data 
lHannam et al\ (|2009b[ ) analyzed some of these gauge dy- 
namics, and developed a promising approach ("trumpet" 
data) to mitigate it (see Sec. IIV.F31 for related gauge is- 
sues). However, for simulations lasting for several orbits 
these modest transient effects are generally negligible. 

Most black-hole binary simulation studies are designed 
to represent the astrophysical population of systems, 
which have circularized before the gravitational radi- 
ation becomes observationally significant. These sim- 
ulations begin with circularly inspiralling initial data 
configurations. Even before reliable numerical relativ- 
ity simulations were possible, considerable attention had 
been given to prescribing initial data for these near- 
circular configurations (see ICookl ( 20001) for a review). 
Within the CTS approach, it is particularly natural to 
impose initially circular motion upon the system, result- 
ing i n an initial data prescrip t ion for quasicircular or- 
bits (ICaudill et all 120061: ICookl l2002t ICook and Pfeifferl . 
I2004D . For either the CTS or Brandt-Briigmann data, 
quasicircular parameters may be chosen by constraints 
on either an effective gravitational potential or total 
energy of the system dCaudill et all . 120061 : ICookl 11994 
iGourgoulhon et all |2002| ). For more realistic inspi- 
ralling traje c tories, the PN approximation may be used 
(|Husa et all l2008bl) . or, for higher accuracy, an it- 
erative procedure in volving short numerical evolutions 
(|Pfeiffer et all 12007ft . 



E. Numerically Friendly Formulations of the Evolution 
Equations 

The Einstein evolution equations, which determine the 
time-development of the initial data, form a set of at least 
six coupled, nonlinear propagation equations. The exact 
formulation of these equations depends on the choice of 
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evolved variables as well as how constraint identities are 
incorporated. Although many formulations are possible, 
not all are equivalent from a numerical point of view. 
The choice of formulation can affect how numerical er- 
rors grow in time, and whether these errors can converge 
to zero as the resolution of the computational grid is in- 
creased. 



1. Hyperbolicity and Well-Posedness 

An ideal formulation is well-posed, meaning that: (i) 
for a given set of initial data, a unique solution exists, 
and (ii) solutions depend continuously on perturbations 
of the initial data. To determine whether a given formu- 
lation has the desirable property of well-posedness, we 
first consider it as a matrix equation: 



d f u = Au + /(u,9ju) 



(17) 



where u is a vector, the components of which represent 
the evolved variables such as components of the metric; 
the linear operator A is the principal part, containing the 
highest-order spatial derivatives; and /(u, Uj) contains 
additional terms which may be nonlinear. Certain condi- 
tions on the eigenvalues and eigenvectors of the operator 
A in Eq. p7|) are used to characterize the equations, typ- 
ically as either "strongly hype rbolic" or " weakly hyper- 
bolic" (ICalabrese et all l2002rj iNagy et all 12004 iReulal . 



120041) . 

If the eigenvalues are real and the eigenvectors form a 
complete set (so any solution can be written as a lin- 
ear combination of them), the equations are strongly 
hyperbolic. Then, assuming adequate boundary condi- 
tions, the equations are well-posed and all solutions are 
bounded by a function that grows (at most) exponen- 
tially at a rate independent of the initial data. Strong 
hyperbolicity is a key ingredient for successful black-hole 
merger simulations. 

On the other hand, if the eigenvalues are real but the 
eigenvectors are not complete, the equations are weakly 
hyperbolic. In this case the equations arc not well- 
posed, and they permit solutions which grow at rates 
that depend on the initial data. Weak hyperbolicity 
implies that small numerical errors in the initial data 
may grow at a rate which depends on the resolution. It 
then becomes difficult to show that the numerical solu- 
tion converges to a well-defined continuum solution, and 
at high resolutions the s imulation may become unstable 
(|Calabrese et all l2002al) . 



2. Harmonic Formulations 

The quest for workable formulations of Einstein's equa- 
tions has proceeded along two parallel lines of develop- 
ment. One originated with consideration of "harmonic 
coordinates", so called because the coordinates satisfy 



the wave equation □ 5 a:'' = 0, where U g is the curved- 
space wave operator. In these coordinates, Einstein's 
equations can be written such that the principal part 
resembles a wave equation in terms of the metric as in 
Eq. (|8]| . In this form, Einstein's equ ations are manifestly 
hyperbolic ( Choquet-Bruhatl . Il962t) . However, the har- 
monic coordinate condition is too restrictive for numer- 
ical purposes, so generalized harmonic coordinates were 
eventually developed by introducing a source term into 



the coordinate condition, i.e. ^gX^ = i?' 1 ( Friedrichl 



119851 : iGarfinkld . l2002h . a suitable choice for which pre- 
serves strong hyperbolicity. The subsequent introduc- 
tion of constraint- damping terms, which tend to drive 
the constraints t oward s zero, further ensured stability 
( Gundlach et al ., 2005). This formulation is manifestly 
second-order in both time and space, and has been im- 
plemented numerically as such ( Pretoriusl 120061 ) , but for 
more efficient numerical integrat ion a first-order- i n-tim e 
formulation was also developed (jLindblom et all 120061 ) , 
and is currently being used by some groups. 



3. ADM-based Formulations 

A second line of development originated with the in- 
vention of the ADM fo rmulation of Einstein's equations 
' 19621) . 



(Arnowitt et all 



This formulation was refined 
bv lYork Jr l (|l979h~ who suggested evolving the specific 
quantity of the extrinsic curvature (Eq. (fTO"|0 . The as- 



sociated evolution equations specify the first-order time- 
derivatives for the three-metric jij and the extrinsic cur- 
vature Kij, as well as for the gauge variables a and /3 l . 
However, the ADM formulatio n is wea kly hyperbolic (see, 
e.g., Chapter 5 of lAlcubierrd ( 20081) ). and attempts at 
stable numerical evolutions with it were not successful. 

There are essentially two ways to modify the ADM 
equations that may affect their hyperbolicity (while keep- 
ing them first-order in time and maintaining their 3+1 
character). One way is to add the constraints, which may 
change the principal part of the equations without chang- 
ing the associated physics. The other is to define new 
independent variables. Simple addition of constraints 
failed to give a strongly hyperbolic formulation however. 
Eventually it was found that strongly hyperbolic ver- 
sions could be constructed from the ADM equations by 
the promotion of certain derived qua ntities to the sta- 
tus of independently evolved variables ( Bona e£~aZl . ll995l : 
iKidder et all 120011 : INagy et q/.l . l2004D . 

Although strong hyperbolicity is important for a sta- 
ble formulation, it is not the only property critical for 
accurate evolution. Of the various strongly hyperbolic 
3+1 formulations, one eventually emerged as more suc- 
cessful than the rest. Its development began with the 
observation that the numerical accuracy of 7, the de- 
terminant of the three-metric, and K, the trace of the 
extrinsic curvature, could best be pr eserved by evolv- 
ing th ese two quantities independently ( Nakamura et all . 
Il987j ). Subsequently, to remove the redundancy in evolv- 
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ing the full three-metric and extrinsic curvature, the con- 
formal three-metric 



and conformal traceless extrinsic curvature 

1 



(18) 



(19) 



were substituted ( Shibata and Nakamural Il995l ). To 
eliminate certain terms with second derivatives (which 
contribute to the principal part of this system), a new 
independently evolved variable 



1" = -djf 



(20) 



was also introduced ( Shibata and Nakamural [l995[ ). Fur- 
ther improvements were made shortly afterwards, such 
as the addition of constraints to eliminate more second- 
derivative terms ( Baumgarte and Shapiro! . Il998ft . The 
resulting formulation, now commonly known as BSSN, 
proved to be robustly stable and accurate for numerical 
purposes. Note that the above choices which comprise 
its form were empirically motivated, by identifying and 
eliminating terms which tended to compromise numeri- 
cal accuracy. Som e analytic justification for its success 
was later found by lAlcubierre et 

Ml (12000ft . who showed 
that the BSSN formulation avoids exponentially grow- 
ing modes and most zero-spe ed modes, whi c h acc umu 



late numerical error. Finally ISarbach et al. I (12002ft and 
iGundlach and Martm-Garci'al J2006|) proved that BSSN 
is also strongly hyperbolic, given an appropriate choice 
of gauge. 

A few no teworthy r efinem ents of the BSSN formu lation 
followed. (|Yo et all l2002ft and iDuez et al\ (|2004l ) sug- 
gested adding specific terms proportional to the "Gamma 
constraint" (r l minus its definition) and the Hamilto- 
nian constraint into the evolution equations. These terms 
have the effect of damping the constraints and thereby 
improving stability. 

Later, for e v olution of punctures specifically, 
ICampanelli et al. | (|2006aft evolved the conformal factor 



in the form x = 7 



-1/3 



which vanishes at the puncture 



(an improvement in accuracy over the previous standard 
of <fi = ln(7)/12, which is singular at the puncture). This 
change of variables, however, introduces the potentially 
singular factor 1/x 2 into the evolution equations, neces- 
sitating an arbitrary restriction on the minimum value of 
X- An alternative variable x' — 7 -1 ^ 6 was subsequently 
suggested because resulting app earances of 1 / x' 2 most! 
cancel with factors of 7 1 / 3 



> earances ot l/x mostly 
dMarronetti et all l2008t 



numerics, especially the stability of the simulation. Im- 
portant considerations include how to deal with the ex- 
treme conditions of black holes such as the physical singu- 
larities, the possible coordinate singularities, the strong- 
field gradients, and the dynamical, surrounding space- 
time. The coordinates must accommodate these features 
in a way that is numerically tractable. 



1. Choosing the Slicing and Shift 

An early consideration was how to avoid the phys- 
ical singularity of a black hole through an appropri- 
ate slicing condition. One such choice, maximal slic- 
ing, required solution of a numerically expensive ellip- 
tic equation (lEstabrook et all Il973t iLichnerowicd . fl944l : 
ISmarr and York Jr.[ 1978ft~ Later, an interest in con- 
structing a hyperbolic evolution system led to a gener- 
alizatio n of harmon i c slici ng known as the Bona-Masso 
family ( Bona et all Il995ft . One particular member of 
this family, called 1+log slicing, was found to also be 
singularity-a voiding, like maximal slicing, but at less nu- 
merical cost ( Anninos et all Il995at iBernsteinl . Il993[) . 

Meanwhile, it was recognized that a shift vec- 
tor was required to counter the large field gra- 
dients, or "slice-str etching", incurred in the pres - 
ence of a black hole ( Alcubierre and Brugmaiml l200lh . 
A hyperbolic shift condition called the "Gamma 
driver" was introduced that fulfilled this requ ire- 



ment (lAlcubierre and Brugmarml . l200ll : lAlcubierre et all 
I2003l . l2001bft~ 



d t r- m B\ 

FB\ 



(21) 
(22) 



iTichv and Marronettil . 120071: Ivan Meterl . 12006ft 



where T l = —djj* depends on a conformal three- 
metric 7ij of the evolving spatial slice, B % is an auxil- 
iary variable, f3 l is the shift, and F is some scalar field. 
T]fi is a damping parameter that fine-tunes the growth 
of the shift, which affects the coordinate size of the 
black-hole horizons, which in turn has bearing on the 
required numerical re solution ( Briigmann et al\ , l2008bt 
iGonzalez et all [2009ft . This shift condition also has the 
desirable tendency to drive the coordinates to quiescence 
in synchrony with the physics, for example after a bi- 
nary merges into a stationary black hole. Clearly if F l 
vanishes, B 1 is damped to zero and /3 l approaches a (typ- 
ically small) constant. 



2. Moving Punctures 



F. Gauge Conditions 

The choice of gauge or coordinate conditions, like the 
choice of formulation, has important consequences on the 



Initially these conditions were intended for use with 
punctures that remained at fixed coordinate positions 
(i.e. comoving coordinates), by choosing F to vanish at 
the punctures. However, this led to large gradients be- 
tween the merging black holes since there the metric had 



14 



to vanish as the physical distance contracted. Patholo- 
gies also arose from the twisting of the coordinates as the 
black holes orbited each other. 

These issues, which resulted in large numerical errors 
and instabilities, were eventually resolved with the break- 
through discovery that slight modifications of the 1+log 
and Gamma-driver conditions allowed arbi trary motion 
of the punctures with robust stability ( Baker et all 
l2006bl: ICampanelli et all l2006al) . In particular, in the 
shift condition, a critical alteration was "unpinning" 
the puncture by no longer requiring the factor F to 
vanish at the puncture but rather to remain constant. 
These modifications were subsequently refined to elim- 
inate zero-speed modes, thus preventing the possibility 
of error build-up ( van Meter et all l2006al ). This was ac- 
complished in part by the addition of advection terms: 



dtB* = dt^-^djT 1 
3 



dtp 1 



'/ ,/>" • ^Ur. (23) 
(24) 



It is noteworthy that these particular gauge choices, 
coupled with the BSS N equations, also result in 
stron g hyperbolicit y (|Gundlach and Martin- Garcial 
120061 ; ISarbach et~al\ . l2002h . A further refinement 
resulted from the observation that, assuming that 
ft = B l = f 1 = initially, Eqs. (J23J) and can be 
integrated to yield B l = Y l — ^r/pft; this substitution 
allows the removal of the equation for B l , lea ving only a 
single shift equation ( van Meter et all l2006aj) : 



dtp 



8 j d j 1' <,hJ 



(25) 



Use of this or similar gauge conditions became known as 
the "moving puncture" method, and proved to be very 
successful as it became increasingly widespread among 
the numerical-relativity community. Analysis has veri- 
fied that moving punctures are valid black-hole solutions, 
although the initial character of the punctures changes 
significantly during the course of evolution. The second 
worldshcct shown in Fig. [31 for example, invariably be- 
comes disconnected in numerical simulations, at a point 
within the horizon near the throat, due to the action 
of the shift vector effectively shifting computational grid 
points onto the first worldsheet. Meanwhile the spatial 
coordinates evolve such that the singularities in the 
conformal factor ip become r -1 / 2 singularities. For a sin- 
gle nonspinning black hole, the numerical result rapidly 
asymptotes to an exact form of the Schwarzschil d solu- 
tion c a lled a "trumpet" , recently investigated by Brownl 
(l2008h:lHannam qZ.I <l2007aL l2009bL l2008bL l2007bT) and 
iBaumgarte and Naculichl ( 2007t ). which turns out to be 
a typ e of solution first considered by lEstabrook et all 
(|!973h . 



3. Generalized Harmonic Coordinates 

Development of generalized harmonic coordinates 
initially proceeded independently of the above 3+1- 



formulated conditions. As mentioned, in harmonic coor- 
dinates the D'Alembertian of each coordinate vanishes. 
In generalized harmonic coordinates, the wave equation 
for each coordinate is allowed a source term, i.e. 



= Hi- 



(26) 



These "gauge driving" source terms can be ei- 
ther algebraically spe cified or e v olved such that hyper- 
bolicity is pres e rved ( Friedrichl Il985t iGarfinklel l2002t 
iLindblom et al. I. I2006t iPretoThisL 120061 ). 

The first successful numerical orbit of black holes in- 
volved a source term for the time coordinate that ef- 
fectively kept the lapse close to its Minkowski value of 
unity, whi l e the spatial coordinates remained harmonic 
(|Pretoriusl . l2006fh This was accomplished by evolving 
the source term itself, according to 



DH 



-Zi(a-l)+&(d t -/3%)H ] a 



(27) 



where £i and £2 are constants. More recently, to 
dampen extraneous gauge dynam ics during the inspi- 
ral and merger of generic binaries, ISzilagvi et all (|2fl09h 
found the following gauge driver to be successful: 



Ho — ^0 
Hi = —/.to 



a 



loe 



JV9 



(28) 
(29) 



where fig is a specified function of time that starts at zero 
and eventually increases monotonically to unity. 



4. Other Coordinate Techniques 

Additional noteworthy advances in coordinate condi- 
tions were developed for the generalized harmonic formu- 
lation to facilitate spectral methods (see Scc. lIV.Gl below) 
but are in principle generalizable to other frameworks. 
One of these is the use of multiple coordinate patches, 
where the coordinates are chosen according to the lo- 
cal physical geometry for optimal accuracy. In the gen- 
eralized harmonic formulation, characteristic speeds are 
readily available for use in constr ucting physical bound- 
ary conditions be t ween patche s (ILindblom et~al. I, l2006t 
iPazos et al. I. I2009t IScheel et al. 1. 120061) . 

The other advance is that of "dual coordinates". Be- 
cause moving punctures involve irregular fields, they are 
not easily made compatible with spectral methods, which 
are sensitive to irregularities. An alternative, that of a 
moving, excised region, can leave a large trail of interpo- 
lation error in its wak e. The remaining optio n is to use 
comoving coordinates (jBriigmann et all 120041 ) . but these 
can result in instabilities due to the steep field gradients 
and other pathologies mentioned in Sec. IIV.F1 Dual co- 
ordinates were invented to exploit the advantages of both 
comoving and non-comoving coordinates while avoiding 
their disadvantages. This is accomplished by computing 
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all tensor field components in the non-comoving coordi- 
nate system (thus avoiding steep gradients), yet evolv- 
ing them as functions of the comoving coor dinates (thus 
allow ing stationary excision boundaries) (jScheel et all 
120061) . 



G. Numerical Approximation Methods 

The initial data, formulation, and gauge are all, in 
principle, analytic choices applicable to an infinite, con- 
tinuous manifold. One must finally make choices pertain- 
ing explicitly to the finite, discrete mechanics of the com- 
putations that will numerically approximate the above 
analytic specifications. An immediate consideration is 
the fact that the simulated domain must have finite ex- 
tent. Various conformal compactification schemes which 
map sp atial or nul l infini t y to a finite boundary have been 
tested <|Pretoriust 120061: iRinnd . l2010t | van Meter et all 
200Gb), but currently it is more common to impose ar- 
tificial boundaries at finite spatial coordinates and ap- 
ply some form of eithe r radiative boundary condition 
(|Alcubierre et all 12003ft or constrain t-preserving bound- 
ary condition (jLindblom et all. 120061) . 



To mitigate the effect of inward-propagating errors 
from the outer boundaries, some form of mesh refine- 
ment is typically employed to push the outer boundaries 
as far away as possible from regions of interest (e.g. wave 
sources and extraction regions). With a limited number 
of grid points available, their density is judiciously chosen 
to be highest near the strong field gradients of the black 
holes, moderate throughout regions where wave propaga- 
tion is studied, and coarser beyond that. If the simulated 
black holes are very dynamic then some algorithm for au- 
tomatically adapting the mesh refinement is necessary. 
Meanwhile the interfaces between refinement regions re- 
quire interpolation. 

On this grid, spatial derivatives are computed in one of 
two ways. They may be computed with finite differencing 
stencils across uniformly spaced grid points, derived from 
Taylor expansions, currently up to eighth-order acc urate 
(|Lousto and Zlochowerl 120081: (Pollnev et all 12009ft . Al- 
ternatively, the spectral approach may be used, in which 
coefficients of an expansion in basis functions are com- 
puted to some order, on a number of collocation points 
comparable to the number of basis modes, fr om which the 
deriva tives are then obtained analytically ( Boyle et all 
l2007bl) . In the former case, dissipative terms are of- 
ten added to the evolut ion equations to reduce noise 
dKreiss and Oligerl I1973D . [n the latter case, spectral 
methods often include a smoothing step for a similar end. 

Last, time must be advanced in discrete steps. Various 
explicit time-integration algorithms have been tested, for 
example the iterated Crank-Nicholson scheme, but the 
fourth-order Runge-Kutta algorithm has become most 
widely used, due to its superior accuracy and efficiency. 
In some codes the timestep size is made to vary with spa- 
tial resolution, for even greater efficiency, albeit at the 



expense of the complication of time-interpolation. Im- 
plicit timestepping schemes are also being inv estigated as 
a me ans to greatly increase the timestep size ( Lau et all , 
I2009D . 



H. Extracting the Physics 

One of the most important end-results of a simulation 
of merging black holes is a computation of the emitted 
gravitational radiation. For this purpose, it is useful to 
calculate the "Weyl tensor", C a bcd- The Weyl tensor, 
like the Einstein tensor, is constructed from derivatives 
of the metric. Unlike the Einstein tensor, the Weyl tensor 
has degrees of freedom that do not necessarily depend on 
a massive source; it can be nonzero while the Einstein 
tensor vanishes. 

Certain components of the Weyl tensor form a com- 
plex quantity called ^4, one of five "Weyl scalar s" used 
to classify spacetimes ( Newman and Penrose! , [T962'*). In a 
special, "transverse-traceless" , spherical coordinate sys- 
tem, ^4 can be expressed as follows: 



r<f>r(f> 



rtird) ' 



(30) 



where the subscripts {f, 9, </>} denote orthonormal tetrad 
components. In a spacetime with gravitational ra- 
diation this quantity typic ally falls off as ~ 1/r 
( Newman and Penrosel Il962[ ). and can be associated 
with outgoing radiation at sp atial infinity (r — > 00) in 
asymptotically flat spacetimes ( Szekeresl Il965[) . In terms 
of the strain introduced in the last section, this can be 
written 



lim (r^i) = lim —r(h + — ih x ) 



(31) 



where, in terms of a metric perturbation h^ v — g^v—Tj^v, 



hy = hai- 



(32) 
(33) 



What makes ^4 a particularly useful measure of the 
radiation is that to first order in /i M „ it is coordinate- 
invariant. 

Although the linearized radiation interpretation for ^4 
is only strictly valid in the limit as r — > 00, in numeri- 
cal simulations it is often extracted on a sphere of large 
but finite radius. This will approximate the expected 
radiative behavior if \h^ v \ << 1; if spheres of multi- 
ple radii are used then it is also possible to extrapolate 
the results to infinity. Typically 3/ 4 is computed on the 
computational grid points from the metric variables and 
then interpolated onto the spherical extraction surface. 
Recently, Cauchy-characteristic extraction methods have 
been applied to all ow direct evaluation of radiation at 
future null infinity ( Reisswig et all l2009h . 

The product of ^4 with a spherical harmonic func- 
tion is then integrated over the sphere, as it is useful 
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to extract specific modes of the field. The structure of 
the gravitational field is such that it is convenient to ex- 
pand in spin- weight -2 spher i cal harmonics _? Y^ m (0, (j>) 
( Newman and Penrose! . H9661 ITeukolskvl [1973) (just as, 
analogously, it proves convenient to use the spin- weighted 
-1 harmonics in expanding the electromagnetic field). 
Thus: 



EE 



VP 



4im -iXlrn 



(34) 



Unlike electrodynamics, there is no dipole moment of the 
gravitational field in conventional general relativity, so 
the dominant contribution is the £ = 2 quadrupole mo- 
ment. 

Other information that is important in analyzing 
spacetimes relates to the properties of the black holes 
themselves. This is commonly taken from apparent hori- 
zons, surfaces from which all light rays must go inward. It 
is useful to compute this during a simulation, for exam- 
ple by finding surfac es on which the e xpansion of light 
rays are minimized ( Thornburd . 120041 ). When a black 
hole is sufficiently isolated, it is possible to define a mass 
and spin magnitude from data evaluated on the hor izon 
(|Campanelli et all I2009L l2006d: iDrever et all 120031) . as 
well as linear momentum ( Krishnan et all 20071 ). These 
quantities can be used to characterize the physical prop- 
erties of the black holes, and to label the spacetimes when 
comparing with PN theory. 



V. BLACK-HOLE MERGER DYNAMICS AND 
WAVEFORMS 

The final merger of a black-hole binary takes place 
in the dynamical, strong-field regime of general relativ- 
ity. For many years, numerical relativists wondered what 
they would find when they probed this arena using black- 
hole merger simulations. What would the merger wave- 
forms look like? How strongly might artifacts in the ini- 
tial data affect the mergers? How might the effects of 
unequal masses and spins change the picture obtained 
by studying the simplest case of nonspinning, equal-mass 
mergers? Would they uncover any unexpected phenom- 
ena? Recent breakthroughs in black-hole merger simula- 
tions are providing important tools for addressing such 
questions. 

In this section we explore the dynamics of black-hole 
mergers and the resulting gravitational waveforms. We 
start with a look at the Lazarus approach, which pro- 
vides hybrid waveforms by combining analytic methods 
with brief numerical simulations. We then begin a more 
comprehensive discussion of black-hole merger physics, 
starting with mergers of nonspinning, equal-mass black 
holes. Next we consider mergers of unequal-mass holes, 
followed by mergers with spin. Throughout this discus- 
sion we aim to provide a historical context while high- 
lighting the key physical results. 



A. First Glimpses of the Merger: The Lazarus Approach 



In the late 1990s, numerical relativity had advanced 
sufficiently to allow brief simulations of two black holes 
in three dimensions. By "brief" we mean here durations 
of ~ 10M, which is a small fraction of the estimated bi- 
nary orbital period of > 100M near the ISCO. While 
most of the community focused on technical develop- 
ments aimed at extending the duration of these simu- 
lations, a small band of col laborators took a different 
approach ( Baker et aUfeOOOl ). 



Their novel idea was to use full numerical-relativity 
simulations to calculate the strong-field approach to 
merger during ~ 10M of simulation time, and then to cal- 
culate t he remaining evolut ion using perturbative tech- 
niques (|Baker et all l2002ari . They be gan with tradi- 



tional puncture black holes (which remain fixed in the 
grid) on quasicircular orbits near the ISCO, and evolved 
them towards merger. Then, just before the simulation 
became unstable, they stopped the calculation and ex- 
tracted the physical data for the merging black holes and 
the emerging gravitational waves. This data was then in- 
terpreted as initial data for a highly distorted single black 
hole, and evolved using techniques of black-hole pertur- 
bation theory. They called this method of reviving a 
(nearly) dead calculation the Lazarus approach. 

The hybrid waveforms produced by the Lazarus col- 
laboration gave the first indications of what might be 
expected from the final m erger of black-hole binaries 
(jBaker et all l200ll . l2002bl ). They ran a suite of simu- 
lations, varying the initial black-hole separations and the 
time at which they stopped the calculation and made 
the transition from the fully numerical evolution to the 
perturbed black hole evolution. Figure 2] shows the dom- 
inant quadrupole £ = 2,m = 2 component of Re(^4), 
which corresponds to the + polarization, for the case 
of equal-mass, nonspinning black holes. Note the sim- 
ple shape of the waveform smoothly tying together what 
might be the end of an (inspiral) chirp with a damped 
(ringdown) sinusoid. 

By tuning the mass and spi n of the background hole 
for the perturbative evolution, iBaker et all (|2002bt ) also 
determined that the merger remnant was a Kerr (spin- 
ning) hole of mass Mxcrr ~ 0.97Mi n i t i a i, and spin axon- ~ 
0.7AW- 

The Lazarus method was also applied to mergers of 
black holes with spins either ali gned or anti-aligne d with 
the orbital angular momentum ( Baker et all 120041 ); in all 
cases, a similar waveform shape was seen. How generic 
was this simple shape and, in particular, would it also 
arise in situations where the black holes complete one or 
more orbits before merging? 
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FIG. 4 Lazarus wa veforms for equal-ma ss, nonspinning black- 
hole mergers, from lBaker et al\ l|2002bt ). The real, I = 2, m = 
2 part of r^A, corresponding to the + polarization (measured 
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each with two transition times to perturbative evolution. 
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B. Mergers of Equal-Mass, Nonspinning Black Holes 

1. The First Merger Waveforms 

Early in 2005, Frans Pretorius stunned the community 
by achieving the first robust simulation of two equal- mass 
black hole s through a sing le plunge orbit, merger, and 
ringdown (|Pretoriusir2005ai) . The resulting gravitational 
waveforms are shown in Fig.[5j where rRe(^E f 4) is plotted 
versus t/Mo- Here r is the coordinate distance from the 
center of the grid to the sphere on which the waveform 
is extracted (see Sec. lIV.Hf . and Mo ~ 0.5M is the mass 
of a single black hole. (Note that the time axes for all 
other waveform plots in this paper are scaled by either 
the total system mass M or the mass of the final, remnant 
black hole Mf). Here curves are plotted for waveforms 
extracted on four spheres with successively larger radii. 
The curves have been shifted in time so that the wave- 
forms overlap. Note that the overall waveform shape is 
simple. The merger yields a single black hole that is spin- 
ning, with spin parameter a/ ~ 0.70M/, where Mf is the 
mass of the final black hole that forms. 

To carry out these mergers, Pretorius ( Pretoriusl 
l2005bL l2006h used techniques that are very different from 
the more traditional approach based on BSSN and punc- 
tures, used by nearly all other numerical relativists. As 
discussed in Sec. IIV.E.21 he used a generalized harmonic 
formulation in which the Einstein equations are written 
with second-order time and space derivatives. The spa- 
tial domain was compactificd, so that the outer boundary 
of each slice was mapped to spatial infinity. Moreover, 
he excised the black holes (each formed from the collapse 
of scalar field "blobs" ) and evolved their motion across 
the grid using adaptive mesh refinement (see Sees. IIV.CI 
and lPvTGl . 

Pretorius had clearly developed a robust method for 
evolving black- hole binary mergers. In the wake of his 



FIG. 5 The first gravitational waveform for black holes evolv- 
ing through a single plunge orbit, merger, and ringdown, 
achieved by iPretorius (2005a). The waves were extracted at 
four radii, and shifted in time to overlap for comparison. Note 
the amplitudes decrease for larger r, due to lower resolution 
in the outer regions 



remarkable success, many numerical relativists began 
studying his methods. However, another surprise was 
soon to emerge. 

In the autumn of 2005 two research groups, one at 
UTB and the other at Goddard, independently devel- 
oped a powerful technique for evolving mergers of punc- 
ture black holes within the traditional BSSN approach 
([Baker et all l2006bl ICampanelli et all l2006al) . As dis- 
cussed in Sec. IIV.F.21 the standard puncture method re- 
quires the black holes to remain fixed on the numerical 
grid. In this new "moving puncture" approach, simple 
but novel gauge conditions allow the punctures to move 
across the grid, producing accurat e and stable merger 
evolutions ( van Meter et all l2006al) . 

Although the UTB and Goddard codes were both 
based on the BSSN approach, they were independently 
written and featured somewhat different implementa- 
tions of moving punctures. In addition, the Goddard 
code used second-order finite differencing (standard in 
the community at the time). They employed a box-in- 
box fixed mesh refinement to produce the high resolution 
in the region around the black holes needed to compute 
the dynamics accurately, while maintaining adequate res- 
olution in more distant regions and a large enough com- 
putational domain to allow accurate extraction of the 
gravitational waves. The UTB code used fourth-order fi- 
nite differencing, and the "fish-eye" coordinate transfor- 
mation to produce higher resolution around the black- 
hole orbits and lower resolution in the wave extraction 
regions. 

Both groups successfully evolved an equal-mass, non- 
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FIG. 6 The trajectories o f the b lack-hole punctures from a 
run by ICampanelli et al\ (|2006aT l . The individual apparent 
horizons are shown at times t = 0, 10M and 18.8A/. The 
solid circles denote the centroids of the apparent horizons ev- 
ery 2.5M during the evolution. The first common horizon, 
marking the time of merger, forms at 18. 8M, just before the 
punctures complete 1/2 orbit. 



spinning binary through the final plunge orbit, merger, 
and ringdown, and extracted the gravitational wave- 
forms. In these simulations, the black holes completed 
~ 1/2 orbit before merger, defined as the time at which 
a single apparent horizon was formed. Figure [6] shows 
the tracks traced by the black-hole punctures, with 
the apparent hori z ons su perimposed, for the UTB run 
( Campanelli et all , l2006al ). These simulations produced 
plunge waveforms with a simple shape similar to that 
found by Pretorius; compare Fig. [5] from Pretorius's sim- 
ulations with Fig. which shows waveforms from the 
Goddard simulation plus a corresponding Lazarus wave- 
form. 



2. Universal Waveform 

In most cases, astrophysical black-hole binaries spiral 
together over many orbits, radiating away any orbital 
eccentricity in the form of gravitational waves. By the 
time the black holes reach the final inspiral, their orbits 
will be quasicircular. All mergers of such equal-mass non- 
spinning binaries should thus produce the same waveform 
signature, subject only to rescaling with the total system 
mass M . 

Over the years, concerns had been raised within the 
relativity community about the effects of deviations from 
astrophysical initial conditions on the waveforms. For ex- 
ample, spurious eccentricity in the orbits could arise from 
starting conditions that did not approximate quasicir- 
cular inspiral accurately. Moreover, the commonly-used 
conformal flatness prescription (see Sec. IIV.D|) for the 



M/24, r=20M 
M/24, r=40M 
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FIG. 7 Gravit ational wave f orms fr om the black-hole merger 
simulations bv lBaker efal] (|2006bf ). The real part of the £ = 
2, m — 2 mode of rVE^ was extracted from the numerical 
simulation on spheres of radii r = 20M, and 40M for the 
medium- and high-resolution runs. The waveforms extracted 
at different radii have been rescaled by r, and shifted in time 
to account for the wave propagation time between the two 
extraction spheres. Curves are shown for medium (M/24 in 
the finest grid) and high (M/32) resolutions. A waveform 
from a Lazarus calculatio n starting with the sa me initial data 
is shown for comparison (B aker et aZ.L [2002bl 'l . The Lazarus 
waveform shown in this figure is scaled differently from those 
in Fig. H 



initial data is different from the conditions experienced 
by an astrophysical binary; these differences would result 
in spurious gravitational radiation being present in the 
initial conditions for the binary simulations. How would 
factors such as these influence the merger waveforms? 

Having developed robust techniques for evolving black- 
hole mergers, numerical relativists eagerly pursued longer 
simulations with the black holes starting at wider sepa- 
rations, completing more orbits, and producing longer 
wavctrains. The UTB group ran a simulation in which 
the black holes completed nearly 1.5 or bits before the for- 
mation of a single apparent horizon ( Campanelli et all , 
l2006bl ). They observed that the resulting waveform was 
quite similar to those produced by their earlier simu- 
lation that started near the plunge (jCampanelli et all , 
l2006al). with a brief oscill atory signal at the beginning. 
( Campanelli et all , l2006b[ ) anticipated "that the plunge 
waveform, when starting from quasicircular orbits, has a 
generic shape that is essentially independent of the initial 
separation of the binary." 

The Goddard group, simultaneously pursuing this 
same goal, produced the first "universal" waveform 
for the merger of e qual-mass, nonspinning black holes 
(Ba ker et M l2006al) . Jsing approximately quasicircu- 
lar initial conditions, they ran a series of four simula- 
tions, starting the black holes at increasingly wider sep- 
arations. In all, the binaries completed ~ 1.8, 2.5, 3.6, 
and 4.2 orbits before forming a common apparent hori- 
zon; the systems emitted just under 4% of their energy 
as gravitational radiation, and the final black holes were 
spinning with (a/M)f ~ 0.69. To compare the results 
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FIG. 8 Puncture tracks from equal-mass, nonspinning black- 
hole merg er simulations starti ng at increasingly larger sepa- 
rations bv lBaker et al\ ((20063). For clarity, the trajectory of 
only one black hole in each run is shown. The tracks lock on 
to a universal trajectory ~ f orbit before merger (denoted by 



of these models, they chose the fiducial time t = to be 
the moment of peak amplitude in the gravitational ra- 
diation; this typically occurs within a few M of merger, 
ft is worth noting how close these results are to those 
obtained from the Lazarus approach to the same system 
(see Sec. Eg} - 

The binary orbital dynamics of these runs is clearly 
seen in Fig. which shows the tracks of the black-hole 
punctures. Here only one puncture track is shown for 
each binary, with the trajectories oriented to superpose 
at the fiducial t = 0. The tracks show the effects of ec- 
centricity in the early stages of each run; plots of the 
separation r versus time show that the initial eccentric- 
ity decreases as the initial separation increases. As the 
holes spiral together into the strong-field regime, the ec- 
centricity diminishes. The puncture tracks lock on to 
a single universal trajectory, independent of the initial 
conditions, through the final orbit, plunge, and merger. 

Figure [5] shows the corresponding universal gravita- 
tional waveform for equal-mass, nonspinning black holes. 
Specifically, it shows the real part of the I = 2, m = 2 
mode of r^i versus time; for this case, the quadrupolc 
(2, ±2) modes strongly dominate all other modes. Here 
the signals produced by each run have been shifted in 
time so that t = marks the moment of peak radiation 
amplitude. Starting from t ~ — 50Mf, the waveforms 
for the final burst of radiation show nearly perfect agree- 
ment, with differences at the level of ~ 1%. For the 
preceding few orbits, the waveforms agree in amplitude 
and phase to ~ 10%, except for the brief initial bursts of 
spurious radiation (see Sec. HV.Dl) . Overall, the merger 




FIG. 9 The universal waveform (£ — 2, m = 2 mode) pro- 
duced by the f our simulations who se puncture tracks are 
shown in Fig. [8] (jBaker et all l2006al ). The signals have been 
shifted in time so that the peak radiation amplitude occurs 
at t = 0. The agreement among the waveforms is excellent, 
with differences ~ 1% for the final burst of radiation starting 
at t ~ -50M/. 



stage lasts ~ 100M, converting ~ 4% of the initial to- 
tal mass M into gravitational-wave energy. The gravita- 
tional wave released during this burst has a luminosity 
L ~ 10 23 L Q , which is greater than the total luminos- 
ity of all the stars in the visible Universe. For stellar 
black-hole binary mergers, this luminosity will last for a 
few milliseconds, while for MBH binaries it will last for 
several minutes. 

Note that this universal waveform has a simple shape. 
The signal starts with a short inspiral chirp, increases 
smoothly in amplitude during the plunge and merger, 
and then transitions to the damped ringdown sinusoid. 
Overall, the frequency increases monotonically, reaching 
a max imum value that st ays constant during the ring- 
down (jBaker et all l2006al) . 

Since black-hole-merger waveforms will be used to help 
find signals in data from gravitational-wave detectors, 
it was important to verify that various groups consis- 
tently achieve the same results. Figure [10] provides a 
comparison of waveforms computed by Pretorius and the 
groups at UTB and God dard for equal-mass , nonspin- 
ning black- hole mergers ( Baker et all l2007b|) ; here the 
real part of r^'s dominant I = 2, m = 2 mode is 
shown, with all three waveforms aligned in time so that 
the moment of peak gravitational radiation amplitude 
occurs at t = 0. All three simulations evolved through 
the last ~ 1.8 — 2.5 orbits before merger; the Goddard 
run is Rl from lBaker et al. I (|2006al) and shown in Figs. [5] 
and Note that all evolutions produce the same overall 
waveform shape. The Goddard and UTB groups evolved 
black holes with zero spin starting on similar quasicir- 
cular initial orbits, and produced nearly identical wave- 
forms. Pretorius used co-rotating initial conditions which 
impart a small spin (a/M)i^ = 0.08 to each black hole. 
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FIG. 10 Comparison of three waveforms from simulations of 
equal-mass black-hole mergers c omputed by Pretori us and 
the UTB and Goddard groups ijBaker et all . l2007bt l. The 
UTB and Goddard simulations used nonspinning black holes, 
whereas Pretorius' initial conditions produced co-rotating 
black holes, each with a small spin (a/M)i,2 = 0.08. The dif- 
ferences between Pretorius' waveform and the others during 
the ringdown, t/Mf > 25, are consistent with his simulation 
producing a slightly faster rotating final black hole due to 
these initial spins. Reproduced by permission of the Institute 
of Physics. 



This difference is consistent with the slightly higher fre- 
quencies seen in Pretorius' results during the ringdown, 
t > 25M/. 

How robust is this universal waveform to larger 
amounts of orbital eccen tricity and gr a vitati onal radia- 
tion in the initial data? iHinder et all (120081) studied a 
series of equal-mass, nonspinning binaries starting from 
roughly the same initial orbital period and having vary- 
ing amounts of initial orbital eccentricity e. For e < 0.4, 
the initial eccentricity is radiated away, and the binaries 
circularize and begin a universal plunge at t ~ 50M/ be- 
fore the time of peak radiation amplitude, producing a 
final black hole with a ~ 0.69A//. For e > 0.5, the black 
holes do not complete any orbits, but rather plunge to- 
gether and merge; the final black hole spin reaches a max- 
imum value a ~ 0.72Mf around e ~ 0.5. Figure ITTI shows 
the puncture tracks, and Fig. [T^] shows the real part of 
the strain's I = 2, m = 2 mode, for two nonspinning, 
equal-mass black- hole mergers with nonzero eccen tricity 
from this study; in both cases. IHinder et al. (2008) show 
quantitatively that the waveform for the final burst of 
radiation is the same as the universal waveform. 



In a complementary study, ISperhake et al. | |2008a) 
found that the final spin of the remnant black hole is 
essentially insensitive to the eccentricity for bi naries that 
do not plunge immediately. iBode et a l. ( 2008) examined 
the effects of spurious gravitational waves by superposing 
a tunable packet of gravitational radiation on an equal- 
mass, nonspinning binary. They found the binary evolu- 
tion and the spin of the remnant black hole to be robust 
to modest amounts of added gravitational radiation. 




FIG. 11 Puncture tracks for equal-mass black-hole mergers 
with initial orbital eccent r icity e = 0.1 (left) and e = 0.3 
(right) from IHinder et al\ l|2008f ). In these cases, the initial 
eccentricity is radiated away and the binaries circularize be- 
fore merger. 
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FIG. 12 Real part of (2, 2) strain mode (labeled h+ here) fo r 
black- hole mergers shown in Fig. llll from lHinder et all (|2008T ). 
These runs have initial orbital eccentricity e = 0.1 (top) and 
e = 0.3 (bottom) and enter a universal plunge by t ~ 50M/ 
before the gravitational radiation reaches its peak value. 



3. Longer Waveforms 

The black-hole simulations discussed so far cover the 
last < four orbits before merger. These achievements 
marked the triumph of solving a long-standing prob- 
lem of fundamental importance to general relativity: the 
two-body problem for the final merger of equal-mass 
Schwarzschild black holes driven by gravitational radi- 
ation reaction. These simulations provided the first look 
at the dynamics and waveforms of black holes evolv- 
ing and merging in the nonlinear, strong-field regime. 
They also pointed the way towards applications in de- 
tecting t he final gravitational- wa ve burst from black-hole 
mergers ( Baumgarte et q/] .l2008l). However, for more ad- 
vanced applications to gravitational-wave data analysis 
and comparisons with PN analytical treatments, longer 
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waveforms that start many cycles before merger are es- 
sential. 

The Goddard group (jBaker et al ], l2007d ldh produced 
the first such long waveforms for equal-mass, nonspinning 
black holes starting ~ 7 orbits or ~ 14 gravitational- 
wave cycles before merger. They used improved initial 
conditions with low eccentricity e < 0.01, and focused 
on improving accuracy while the black holes traversed 
the relatively long inspiral. They investigated the ob- 
servability of black-hole mergers wit h ground- and space - 
based gravitational-wave detectors (|Baker et all l2007cf ) 
(see Sec lVII.A| and successfully applied their long wave- 
forms to compari sons with PN result s, focusing on the 
waveform phases ( Baker et all l2007dl ) . 4 

Shortly thereafter, the Jena group simulated a bi- 
nary inspiralling for nine orbits (18 gravit ational- wave 
cycles) before merger ( Hannam et oZ.I . [2008ch . With this, 
they made the first quantitative comparisons with both 
the PN phase and amplitude, and quantified the level 
of error in the quadrupole app roximation. They used 
higher-order finite differencing ( Husa et all l2008af ) and 
initial binary parameters calculated using the PN ap- 
proximation to redu ce the initial eccentricity significantly 
( Husa et al 1 l2008bl) . enabling a precise measurement of 
the waveform phase. 

The Caltech-Cornell group currently holds the record 
for the longest and most accurate black-hole binary evo- 
lution, starting 16 orbits and 32 gr avitational- wave cycles 
before merger (jScheel et all I2009T ) . Using their spectral 
code (see Sec. lIV.G|) . they begin with a very small initial 
orbital eccentricity e ~ 5x 10~ 5 and evolve with very high 
accuracy through a relatively long inspiral, then merger 
and ringdown. The impressive trajectories of their black 
holes are shown in Fig. fTS] and the accompanying gravi- 
tational waveforms in Fig. 1141 

Comparison of waveforms from different groups re- 
mains important in the push for higher accuracy and 
use in gravitational-wave data analysis. The Samurai 
project I Hannam et all [2009al) sets the current state of 
the art for studying the consistency of black-hole binary 
waveforms. This effort focuses on comparing waveforms 
from equal-mass, nonspinning binaries, starting with at 
least six orbits (or 12 gravitational- wave cycles) before 
merger and continuing through the ringdown. They fo- 
cus on comparing the amplitude A(t) and phase <f>(t) for 
the i = 2, m = 2 mode of r^, defined as 



r* 4 , 22 (t) = A(t)e-^«; 



(35) 



the gravitational-wave frequency of this mode is then 



4 Recall that the PN approximation is an expansion in powers 
e = v 2 /c 2 and applies when the black holes are far enough 
apart that the black hole speeds remain well below the speed 
of light. We refer to PN results by the order at which the series 
is truncated. For example, "2 PN" means that terms of order 
e 2 = t> 4 /c 4 are retained. See Sec. I VII for a deeper discussion of 
the PN approximation in the context of numerical relativity. 




FIG. 13 Trajectories for the merger of equal-mass, nonspin- 
ning black holes computed by t he Caltech-Cornell g roup using 
their spectral evolution code (|Scheel et all 120091 ). The red 
and black circles/ellipses are the initial/final coordinate loca- 
tions of the apparent-horizon surfaces, and the blue thick con- 
tour is the common apparent horizon just after it appeared. 
The black holes complete 16 orbits before merging. Figure 
kindly provided by H. Pfeiffer. 
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FIG. 14 Gravitational waveforms from the Caltech-Cornell 
merger simulation seen i n Fig. 1 131 showing the £ = 2, m = 2 
component of Re(r* 4 ) (|Scheel et all 120091 ). The left panel 
shows a zoom of the inspiral waveform, and the right panel a 
zoom of the merger and ringdown. 




uj(t) = 4>{i). They compare the results from five in- 
dependent numerical codes: the moving-puncture codes 
from the AEI, Goddard, Jena, and Penn State groups, 
and the Caltech-Cornell spectral code. Figure [15] com- 
pares the gravitational- wave amplitudes and Fig. 1161 the 
gravitational-wave phases as a function of frequency for 
the five waveforms. Qualitatively, the results appear to 
be quite consistent. Quantitatively, they concluded that 
these waveforms have sufficient accuracy to be used for 
detection with all current and planned ground-based de- 
tectors. 
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FIG. 15 The gravitational-wave phase <f> as a function of the 
dimensionless fr equency Mo> for the five codes from the Samu- 
rai comparison ijHannam et all l2009al ) . The scale along the 
top of the panel labels the frequency in kHz scaled with re- 
spect to the total binary mass in solar units. 
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FIG. 16 Amplitudes of the £ — 2, m = 2 component of the 
gravitational waves produced by the five codes in the Samurai 
comparison (jHannam et al ., 2009a) are shown as a function 
of frequency. 



C. Mergers of Unequal-Mass, Nonspinning Black Holes 

Early in 2006, numerical relativists took the next step 
in opening up the parameter space of binary black-hole 
mergers by simulating nonspinning binaries with unequal 
masses. This added a new parameter, the mass ratio q, to 
the problem and brought an additional need for adaptive 
mesh refinement to achieve adequate resolution around 
the smaller black hole. As shown in Fig. [TTJ the smaller 
black hole also moves faster, completing an orbit around 
the center of mass in the same time as the larger hole and 
thus requiring smaller timesteps for its evolution. These 
factors combine to make simulations of unequal-mass bi- 
naries more technically challenging; currently, numerical 
relativists are able to simulate mass ratios up to q = 10 
dGonzalez et all 120091) . 
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FIG. 17 Puncture tracks for an unequal-mass nonspinning 
binary with q = 2 produced by the Caltech-Cornell spec- 
tral code; figure kindly provided by H. Pfeiffer. The solid 
(red) curve gives the trajectory of the smaller black hole, and 
the dotted (blue) curve the path of the larger one. The red 
and blue contours are the final coordinate locations of the 
apparent-horizon surfaces, and the black contour is the com- 
mon apparent horizon just after it appeared. 



1. Mode Analysis and Gravitational Waveforms 

Decomposition into spin-weighted spherical harmonic 
modes (see Sec. HV.Hp provides the basis for an in-depth 
study of black-hole mergers. For the equal-mass case 
5 = 1, the £ = 2, m = ±2 quadrupole mode is dominant 
and the odd-m modes are suppressed by symmetry As q 
increases, the sub-dominant modes become more impor- 
tant, affecting both the evolution of the source and the 
emitted radiation. 

iBerti et all (|2007bt) carried out a set of unequal-mass 
mergers with mass ratios in the range 1 < q < 4. They 
found that, to leading order, the total energy emitted 
during merger scales ~ rj 2 while the spin of the final black 
hole scales ~ 77, where r] is the symmetric mass ratio ([1}. 
They also studied the multipolar structure of the gravi- 
tational waves. Figure H51 shows several (£,m) modes of 
the radiation produced in their simulations for the case 
q = 2. They demonstrated that the higher modes carry 
larger fractions of the total energy as q increases; in par- 
ticular, the £ = 3 mode generally carries ~ 10% of the 
emitted energy for q > 2. 

The Goddard group (jBaker et al I l2008bl ) performed 
a complementary study of nonspinning unequal-mass 
mergers for mass ratios 1 < q < 6. They found that the 
overall simple waveform shape first discovered for equal- 
mass mergers extends to the cases with q > 1; this is 
easily seen when the gravitational waves are scaled by 77. 
Figure [151 shows the strain rh+/r), including the 1 = 2 
and £ = 3 modes, for an observer located at distance r 
on the system's orbital axis. The waveforms are aligned 
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FIG. 18 The amplitudes of several (I, ro) modes of Mr^ 4 
for mass ratio q = 2 from simulations bv lBerti et all (|2007bT) . 
The initial burst of radiation is an artifact of the initial data, 
and the wiggles at late times result from numerical noise. 
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FIG. 19 Strain waveforms sca led by 77 for different mass ra- 
tios, from iBaker et all l|2008bh . h+ is calculated for different 
mass ratios using the £ — 2 and £ — 3 modes. The observer is 
located at distance r along the axis of the system. 



so that the maximum amplitude of the I = 2, m = 2 
mode occurs at t = 0. The differences in phase during 
the final ringdown portion of the waveforms are consis- 
tent with the differing spins of the final black holes; for 
example, with nonspinning black holes, a larger mass ra- 
tio q results in a smaller final spin for the remnant black 
hole. When viewed off-axis, the waveforms show modest 
amplitude and phase modulations, while still preserving 
the overall simplicity in shape. 

IBaker et al. (2008b) also found that each of the indi- 
vidual spherical harmonic components is circularly polar- 
ized during the inspiral, merger, and ringdown, as seen 
by distant observers on the system's rotational axis. Dur- 
ing the inspiral, the phase and frequency of the different 
(£, m) components are the same for each mass ratio; this 
is expected since the waveform phase is directly related 
to the orbital phase. More interestingly, for the I = m 
modes, this relationship continues to hold during merger 
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FIG. 20 The £ = m = 2a,nd£ = m = 3 compo- 
nent of Re(rMf 4) fr om the 10:1 mass ratio simulation by 
iGonzalez etoH (20(jl) 



and ringdown. 

Based on these observations, they developed a simple 
conceptual picture in which each (£ , m) mode of the grav- 
itational radiation is produced separately by the (£, m) 
mode of some implicit rotating source. The fixed rela- 
tionship between the phase and frequency of the £ = m 
modes shows that these components of the implicit source 
rotate rigidly during the entire coalescence from inspiral 
through ringdown. In comparison, the I 7^ m source 
components peel away from the main trend of the £ = m 
p arts during merger, ind icating less rigid rotation. 

I Gonzalez et ~al. (|2009h have carried out the highest 
mass ratio merger simulation to date, with q = 10. This 
binary radiates <~ 0.42% of its mass as gravitational ra- 
diation as it undergoes ~ 3 orbits before merging to 
form a black hole with (a/M)f ~ 0.26; th e se valu es fit 
the scalin g relat ions found by iBerti et~al\ (l2007h|) and 
iPan et all (|2008l ). Figure!^ from lGonzalez et all (|2009h 
shows that the £ = m = 2 and I = m = 3 modes of 
Re(rM$4) exhibit the simple waveform shape individu- 
ally; this trend continues through £ — m = 5, the highest 
mode they studied. 



2. A Qualitatively New Feature: Kicks 

Unequal-mass black-hole binaries can form as a result 
of galaxy mergers, dynamical processes in star clusters, 
and the final stages of binary stellar evolution. Mergers 
of unequal-mass binaries bring a qualitatively new phe- 
nomenon of great importance to astrophysical scenarios 
of black-hole growth and retention: recoils or kicks. 

As discussed in Sec. III.Al the gravitational waves emit- 
ted by a merging black-hole binary carry away linear mo- 
mentum. If the pattern of gravitational- wave emission is 
not symmetrical {i.e., if there is more radiation emitted 
in some direction than in others), then global conserva- 
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FIG. 21 Schematic drawing of the physical basis of the recoil 
or kick produced by a merging black-hole binary with unequal 
masses. A net flux of momentum P e j ec ted is emitted parallel to 
the smaller hole's velocity. Momentum conservation requires 
that the system center of mass recoil in the opposite direc- 
tion, Prccoii- Reproduced with kin d permission fr o m Sp ringer 
Science+Business Media: Fig. 1 of iHughes et all (|2005h . 
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FIG. 22 The magnitude of the kick velocity fo r a nonspin- 
ning, u nequal mass merger with q — 1.5 from iBaker et all 
(|2006d ). Results are shown from three simulations, with 
increasingly wider initial black-hole coordinate separations 
di n it. The small-dotted (purple) curve shows a 2 PN calcu- 
lation starting with initial conditions commensurate with the 
dinit = 7.0M case (blue dotted line). Reproduced by permis- 
sion of the A AS. 



tion of momentum requires the center of mass, and thus 
the remnant black hole that forms, to recoil in the oppo- 
site direction. 

The situation for an unequal-mass, nons pinning binary 
i s sho wn schematically in Fig. from IHughes et all 
(I2005D . who attributed this argument to Alan Wiseman; 
we summarize their discussion here. The two holes are 
orbiting in a plane about their center of mass. Since the 
smaller black hole moves faster, its wave emission un- 
dergoes more "forward beaming" than that of the larger 
hole. Instantaneously this gives a net flux of momentum 
parallel to the smaller hole's velocity, and an opposing 
recoil or kick at the center of mass. The direction of this 
kick changes continually as the black holes orbit. If the 
orbit was circular, the center of mass would move in a 
circle and suffer no net recoil. However, since the black 
holes are spiralling together (due to the energy and angu- 
lar momentum carried away by the gravitational waves) , 
the instantaneous kicks do not cancel exactly but rather 
accumulate. This causes the final merged black hole to 
have a non-zero net recoil in the orbital plane. 

This recoil has been studied by several au- 



thors, including Bekensteinl (11973 ); iFitchettl 



(119831); iFitehett and Detweilerl dl984l): iPeresI (Il962l) 



iRedmount and Reesl ( 19891 ). More recent analytic 
treatments using PN a pproximations were carried out b y 



treatments using rlN approximations wer e carried out b y 
Blanchet et dl (120051) : I Damour and Gopakumarl (120061); 



and ringdown ( Schnittman et all , 120081 ) . 

The Goddard group (|Baker et all , l2006ch carried out 
the first accurate calculation of black-hole recoil for the 
merger of a nonspinning binary with q = 1.5. Figure [2"2l 
shows the kick velocity as a function of time for three 
simulations in which the black holes start out on orbits 
with increasingly larger initial separations d; n it. The kick 
velocity grows rapidly during merger, reaching a peak 
value ~ lOM after the time of peak radiation amplitude 
in \l/4, and then dropping to a lower, final value. It is 
clear from this figure that the black holes must start far 
enough apart in order to get a consistent and reliable 
value for the final recoil velocity. For this q = 1.5 case, 
their results spanned the range uidck = 101 ± 15 km s~ , 
with a best estimate of fkick = 92 ± 6 km s _ . 

Following on from this, and other u n equal-mass sim- 
ulation s by iHerrmann et ~aH (|2007bl ). IGonzalez et~al\ 
( 2007bl ) performed the first systematic parameter study 
of nonspinning black-hole binary mergers with mass ra- 
tios in the range 0.253 < 1/q < 1. They found that the 
maximum recoil velocity Ukick = 175.2 ± 11km s _1 oc- 
curs for the mass ratio 1/q = .36 ± 0.03. For the case 
q = 10, IGonzalez et al. I (12009D found a recoil velocity 
^kick = 66.7 ± 3.3 km s _1 . 



(|2004jk iLe Tiec et all (|2010D ; IWisemanl D. Mergers of Spinning Black Holes 



Favata et al. 

(11992D . while some comparable-mass est imates were 
also made with the Lazarus approach (|Campanellil . 
I2005D . However, since the dominant part of the kick 
depends sensitively on the strong-field regime close to 
merger, where the orbits are less circular, an accurate 
calculation of the kick requires full numerical relativity 
simulations of the final stages of binary inspiral, merger 



Astrophysical black holes are generally expected to be 
spinning. Including the spin vector of both black holes 
in the binary adds six more dimensions to the param- 
eter space, giving seven in all when the mass ratio q 
is included. Exploration of this large parameter space 
began in 2006, with the simplest cases of equal masses 
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FIG. 23 Gravitational waveforms for mergers o f equal-mass 
spinning binaries as in ICampanelli et al\ (|2006dl ) . The domi- 
nant I — 2, m = 2 mode is shown for three cases: aligned spins 
(solid/red), anti-aligned spins (dotted/black), and no spins 
(dashed/blue). All three simulations start with the same ini- 
tial orbital period at t = 0. Figure kindly provided by M. 
Campanelli. 



and spins, and is gradually growing to incorporate more 
generic binaries. 



1. Gravitational Waveforms 

The UTB group carried out the fir s t black -hole binary 
mergers with spin ( Campanelli et alX l2006dl) . They sim- 
ulated three binaries, each having equal masses and start- 
ing on quasicircular orbits with initial orbital frequency 
MO = 0.5, giving an initial orbital period ~ 125M. 
In the aligned case, both holes have spin parameters 
(a/M)i,2 = 0.757 and vectors parallel to the orbital an- 
gular momentum L; in the anti-aligned case, the black 
holes have the same spin parameters but the vectors are 
anti-parallel to L. A nonspinning equal-mass binary was 
run for comparison. 

The I = 2, m = 2 components of the waveforms are 
plotted in Fig. [531 which shows all three cases starting 
out at the same time t = at the same orbital period 
and gravitational-wave frequency ,/gw ~ 2/ or b- Notice 
that the aligned system takes longer to merge, completing 
more orbits and producing a longer wavetrain. This be- 
havior is caused by a spin-orbit interaction that produces 
an effective repulsive force between the black holes. In 
the anti-aligned case, this interaction yields an effective 
attraction, causing the binary to merge more quickly and 
with fewer orbits. The nonspinning case is intermediate 
between the other two. All three cases produce similar 
waveforms with a clean, simple shape. 

It is also interesting to compare the total energy radi- 
ated as gravitational waves £" ra d and the spin parameter 



of the final black hole (a/M) / for these three cases. The 
aligned case yields E lad ~ 0.07M and (a/M)/ ~ 0.89 
while the anti-aligned case gives -Erad ~ 0.02M and 
(a/M)/ ~ 0.44. The nonspinning case is again between 
these two cases, w i th E ra ,d ~ 0.04M and (a/M)/ ~ 0.69 
(jCampanelli et aU |2006cl . 

The black-hole mergers discussed so far have resulted 
from simple binary dynamics with no precession and have 
produced simple waveforms with a similar shape, at least 
when examined mode-by-mode in a spin- weighted spher- 
ical harmonic decomposition. Might more "generic" pro- 
cessing binaries generate waveforms with more complex 
patterns? This important question pertains not only to 
basic orbital dynamics in general relativity, but also to 
strategies for developing templates to search for signals 
in data from gravitational- wave detectors (see Sec. IVII|) . 

The parameter space for such binaries is very large, and 
explorations of this space have only recently started. The 
UTB group, now at the Rochester Institute of Technol- 
ogy (RIT), pioneered the study of precessing black- hole 
binary mergers. They evolved a generic black-hole binary 
with unequal masses and unequal, non-aligned precess- 
ing spins that undergoes ~ 9 orbits before merger and 
prod uces a relatively long wavetrain ( Campanelli et all , 
12009t h This system has a mass ratio q = 1.25 and arbi- 
trarily oriented spins with magnitudes (a/M)i ~ 0.6 and 
(a/M)2 ~ 0.4. The initial conditions were determined 
from point-mass evolutions using 2.5 PN and 3.5 PN pa- 
rameters. 

Figure [2U shows the difference in the black-hole tra- 
jectories x\ — X2- For a non-precessing binary, x\ — x*i 
would have no component in the z direction. In this 
generic case, the precession of the total system spin drives 
a precession of the orbital plane, producing evolution of 
the trajectory in the z direction. The resulting wave- 
forms demonstrate the effects of precession on the am- 
plitudes of the sub-dominant modes. Figure [55] shows 
the I — 2, to = 1 mode of the strain h+. The numerical- 
relativity evolution is shown by the solid curve; note the 
amplitude modulations induced by the precession. The 
dotted curve shows a solution of the PN equations of 
motion at 3.5 PN order. 



2. Spinning Binary Mergers and Spin Flips 

The merger of a binary consisting of nonspinning black 
holes produces a single spinning black hole with spin di- 
rection aligned with the orbital angular momentum L 
of the binary. In this case, the spin of the final hole 
arises purely from L. When the individual black holes 
have spins that are not aligned with L, the spin of the 
final hole will, in general, not be aligned with the initial 
spin but rather will have a different direction. This phe- 
nomenon, in which the spin direction of the final black 
hole differs significantly from the spin directions of the 
individual holes prior to merger, is known as a "spin flip." 
The simplest situation in which a spin flip can occur is a 
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FIG. 26 Black hole trajectorie s for the spin flip case stud- 
ied by ICampanelli et al\ (|2007al ). The black holes have equal 
masses, and equal spins initially pointing 45° above the initial 
orbital plane. The spin directions are shown as arrows along 
each trajectory every 4M until merger. 



FIG. 24 The difference in the black -hole trajectories x\ — x *2 
for the generic binary evolution from lCampanelli et al\ (|2009t ). 
The precession of the system spin induces precession of the 
orbital plane. 
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FIG. 25 Precession-induced amplitude modulations in the 
waveforms from the gen eric black-hole binary merger from 
ICampanelli et all (|2009h . The 1 = 2, m = 1 mode of the 
strain (h+) is shown for the numerical-relativity simulation 
(solid) and a 3.5 PN calculation (dotted). The amplitude 
oscillations are induced by precession. 



binary with both spins anti-aligned with L. 

The case of generic binaries, with unequal masses and 
arbitrarily-oriented spin s has a much larger par ameter 
space. The RIT group ( Campanelli et a/.U2007al lbl) has 
examined spin flips for several such cases. Figure [26] 
shows the puncture trajectories for an equal-mass binary 
with equal spins having magnitude (a/M)i )2 = 0.5007 
and initially pointing 45° above the initial orbital plane. 
The black holes execute ~ 2.25 orbits before merger, dur- 
ing which time the spins, shown as arrows along each 
trajectory, precess by ~ 151°. The evolution of the spin 



direction for this case is shown in Fig. [27] The final black 
hole has (a/M)f ~ 0.8 with a spin direction flipped by 
~ 35° with respect to the component spins just before 
merger. 

In a related phenomenon, the direction of the total 
angular momentum [L + S\ + Si)) may change. This 
case was studied using general princi ples and extrap- 
olatio ns from point-particle results by iBuonanno et gj\ 
( 20081 ). who predicted that binaries with q > 6.78 will 



experience a flip in total angular momentum direction 
provided that the initial spins are equal and anti-aligned 
with L, and that the individual black hole spins obey 
(a/M)i t 2 > 0.5. In particular, they predicted that, 
for q = 4, a Schwarzschild black hole should form if 
(a/M) 1 = {a/M) 2 = 0.8. 

Motivated by this work. lBerti et al. I (I2008D simulated a 
series of binaries with q — 4 and having equal anti-aligned 
spins in the range (a/M)i, 2 € [-0.75, -0.87]. They found 
that a Schwarzschild hole is formed when the initial spin 
is {a/M)i, 2 ~ -0.842 ±0.003. 



3. Kicks from Mergers of Spinning Black Holes 

The mergers of asymmetrical spinning binaries will 
produce recoiling remnant black holes. In the simplest 
case, the binary components have equal masses and spin 
magnitudes, and the asymmetry arises from the spin vec- 
tors pointing in opposite directions. In early 2007, several 
new results appeared in quick succession, generating con- 
siderable excitement among numerical relativists as well 
as astrophysicists eager to apply these results. In this 
section we discuss the basic results and defer to Sec lVIHl 
t heir impact on as tr ophysic s. 

iHerrmann et al. I (l2007aD simulated the merger of 
equal-mass black holes with equal spin magnitudes 
(a/M)i i2 € (0.2,0.4,0.6,0.8) and having one spin vec- 
tor aligned and the other anti-aligned with the orbital 
angular momentum. These mergers produced kicks di- 
rected purely in the orbital plane with magnitudes ~ 
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FIG. 27 The spin flip produce d by the evolution in Fig. [26] 
from ICampanelli et all lj2007aT ). The spin direction of the 
component holes, shown by the red arrows plotted every AM 
until merger, varies continuously due to precession. The blue 
arrow shows the final black-hole spin, with direction flipped 
discontinuously from that of the component holes just prior 
to merger. 



475(a /M)i .2 km s 1 . Shortly thereafter iKoppitz et all 
( 20071 ) showed that, for the specific case (a/M)i^ = 
0.584, Vkick = 257±15km s _1 ; linear extrapolation of re- 
sults to the maximally spinning case (a/M) 1: 2 = 1 yields 
w kick ~ 440 km s _1 . iPollnev et all (|2007l ) carried out a 
systematic study of equal-mass mergers with spin mag- 
nitude (a/M) 2 = 0.584 and direction aligned with the or- 
bital angular momentum. The other black hole has spin 
magnitude (o/M)i G (0, 0.25, 0.50, 0.75, 1.0)(a/M) 2 , and 
direction both aligned and anti-aligned. They found a 
maximum kick velocity Wkick = 448 ± 5 km s" 1 . 

These results all demonstrated that mergers of spin- 
ning black holes can produce significantly larger kick ve- 
locities than nonspinning mergers. However, new results 
would soon show that kicks from the mergers of spinning 
black holes could get much larger indeed. 

The idea of "superkicks" was fi rst discussed by the RIT 
group ( Campanelli et al . 2007bl ). They observed that a 
PN treatment ( Kiddeii 19951) predicts that the recoil due 
to spin is maximized when Mi = M 2 , (a/M)i = (a/M)^, 
and the spin directions are anti-aligned with each other 
and lying in the orbital plane. 

The Jena group first simulated binaries in this con- 
figuration and demonstrated the phenomenon of super- 
kicks. Using (a/M)i i2 ~ 0.8, they found a resulting kick 
velocity Vkick ~ 2500km s _1 ( Gonzalez et all l2007a[) . 
Figure [35] shows the coordinate positions of the black- 
hole punctures from one of their simulations. Notice 
that the trajectories move out of the initial orbital plane 
and that the final bla ck hole is kicked in the — z direc- 
tion. The RIT group ( Campanelli et aZl . |2007cT ) carried 
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FIG. 28 Punc ture trajectorie s from a superkick simulation 
carried out bv lGonzalez et all (|2007al ). The black holes have 
equal masses and spins (a/M)i,2 ~ 0.8, initially oppositely 
directed in the orbital plane. During the evolution, the black 
holes move out of the original plane, and the final black hole 
recoils with velocity «kick ~ 2500 km s _1 in the negative z- 
direction. 



out similar simulations with (a/M)i£ = 0.5. Their re- 
sults show kicks perpendicular to the orbital plane with 
magnitud es up to y k j c k ~ 1800 km s _1 ; using the expres- 
sion from iKidderl ( 19951 ) they predict a maximum recoil 
fkick ~ 4000 km s" 1 for the case of maximally spinning 
bl ack holes, (a/M) 1 , 2 = 1. 

ISchnittman et al\ (|2008f ) performed a multipolar anal- 
ysis of recoil from black-hole mergers, for both unequal 
masses and non-zero, non-precessing spins. They found 
that multipolc moments up to and including £ = 4 
are sufficient to accurately reproduce the final recoil 
velocity (within ~ 2%), and that only a few domi- 
nant modes cont r ibute s ignificantly to it (within ~ 5%). 
iBriigmann et all (|2008al ) studied the role of spin in pro- 
ducing superkicks. They showed that the recoil veloc- 
ity is almost entirely due to the asymmetry between the 
(t = 2,m = +2) and (£ = 2,m = -2) modes of * 4 . 
In addition, the major contribution to the recoil occurs 
within a period ~ 30Af before and after merger, after the 
time at which a standard PN treatment of the evolution 
breaks down. 



VI. INTERACTION OF NUMERICAL RELATIVITY WITH 
POST-NEWTONIAN THEORY 

Einstein's theory of general relativity implies predic- 
tions for the dynamics and gravitational-wave genera- 
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tion for generic black-hole binaries at all stages of the 
merger process. As we have seen, numerical relativity 
now provides an excellent tool for concretely deriving 
these predictions beginning in the last orbits and con- 
tinuing through the conclusion of the coalescence in the 
ringdown of the remnant black hole. There is no hard 
limit to how early numerical relativity can be applied 
in these simulations, but the considerable costs of these 
simulations tends to limit the duration to tens of orbital 
periods. Though the PN expansion provides only an ap- 
proximation to general relativity, it succeeds in making 
accurate and efficient predictions over very long time- 
scales, while the black holes remain well separated. Late 
in the merger, however, as the black-hole velocities in- 
crease, these errors grow, and PN theory is no longer 
reliable 

A full understanding of general relativity's predic- 
tions for the merger process must draw on both ap- 
proaches. The maturation of numerical relativity has 
fueled efforts toward synthesizing a complete understand- 
ing of black-hole mergers built from the results and tech- 
niques of both these fields. Even a coarse general re- 
view of the PN formalism and results is outside the 
scope of this presentation. De cades of work in this 
area have already been covered dBlanchet j_ 2006 L 201C ; 
Blanchet efail 120081: iDamour and Nagaiil2010HSchafer , 
2010b . Here we focus exclusively on key synergistic areas 
where fruitful research interactions drawing from both 
numerical relativity and PN results and formalisms are 
yielding a more complete general-relativistic understand- 
ing of black-hole binary systems. The broad range of 
research that interfaces these two approaches can be 
grouped into four categories, based on the different ways 
that PN theory impacts numerical-relativity research. 
PN theory: 1) provides independent results for cross- 
checking and comparing with NR; 2) provides models for 
the initial values needed to begin astrophysically realistic 
NR simulations; 3) provides insight for interpreting NR 
results; and 4) may provide the basis for empirical mod- 
els representing the combined knowledge drawn from PN 
and NR investigations. 



A. Independent Post-Newtonian Dynamics and Waveforms 

As mentioned, PN theory is based on an approximate 
expansion of Einstein's equations in powers of velocity 
e n ~ (v/c) 2n providing general relativistic corrections 
to the Newtonian, small-velocity limit. This approach 
has a long history of success as the primary framework 
for deriving explicit predictions from general relativity 
for physical and astrophysical observations. The rela- 
tively small velocities involved in most observations have 
made low-order PN theory an ideal tool for testing gen- 
eral relativity as the standard model for gravitational 
physics in many contexts including Earth-orbit experi- 
ments, solar-system dynamical observations, and in even 
precision pulsar timing observations of binary neutron 



star systems (| Willi . 120061 ). 

Gravitational- wave researchers have applied PN theory 
to represent general relativity-based signal expectations 
for the vast majority of gravitational- wave searches for 
anticipated observations of black-hole and neutron-star 
binary inspirals, and likewise in the development and 
planning for current and future gravitational-wave in- 
struments. The early epoch of black-hole binary inspiral 
is also well-described by low-order PN theory. For near- 
circular inspirals, e ~ (v/c) 2 ~ M / R, where the binary 
separation R is scaled by the total system mass M. At 
large R, the PN expansion provides excellent predictions, 
even at low order. As the black holes lose energy and 
sink closer together, the velocity grows, requiring higher- 
order PN terms for sufficiently accurate predictions. Cur- 
rently, PN predictions are available up to 3.5 PN order 
(2.5 PN for spin-terms). These higher-order expansions 
seem to provide PN predictions sufficiently accurate for 
the analysis of data from current instruments for all but 
the last several orbits. Once the separation approaches, 
say, R/M ~ 10, the accuracy of the PN expansion di- 
minishes. Even at (hypothetical) arbitrarily high orders, 
the expansion may fail, if the expansion parameter v/c 
exceeds the series' radius of convergence. Precisely when 
the PN approximation effectively fails depends on the 
details of the system being studied and the requisite ac- 
curacy, but generally, for the last orbits and merger PN 
theory is no longer reliable. 

The strongest gravitational radiation is generated in 
the late stages of inspiral or merger where the internal 
consistency of t he PN approa c h has been, at best, dif- 
ficult to assess ( Simone et all Il997t ). Numerical simu- 
lation can treat the late portions of the mergers, with 
practical resource limitations on the duration of the sim- 
ulations, and consequently how far apart the black holes 
can be at the start of the simulation. Would it be practi- 
cal, however, to run numerical simulations long enough to 
"overlap" with the part of the problem treated success- 
fully by PN methods, or would there be an intermedi- 
ate r egion requiring yet another approach ( Brady et all 
Il998h ? 

The first numerical inspiral r esults were quickly 
compared again s t PN calculations ( Baker et all l2006at 
iBuonanno et al. l. l2007aT) . These first comparisons yielded 
promising indications that the gap between NR and PN 
could be bridged, but also made clear that PN results 
were not uniquely determined for the final part of the 
inspiral. Numerical results can be verified by internal 
consistency studies, examining, for instance, the conver- 
gence of the results toward consistency with Einstein's 
equations as the resolution is increased. External veri- 
fication, however, would strengthen the case that these 
new results are indeed correct. Even more importantly, 
these comparisons would allow an independent check of 
the late inspiral PN predictions which are not easily con- 
firmed by self-consistency studies. 

Roughly one year after the first robust numerical- 
relativity results, simulations lasting ~ 1000M were con- 
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FIG. 30 Phase differences between a numerical simulation 



and various PN models from iBovle et all (|2007al ). The left 
FIG. 29 Compari son of NR and PN waveforms from Panel shows phase differences from an initial (waveform) an- 
iBaker et all (|2007dh provided mutual affirmation of the re- gular frequency of Mu = 0.04 up to Mu = 0.063, while the 
suits from each approach and showed that in combination comparison in right panel extends this range to Muj = 0.10. 
NR and PN results could treat the complete signal. 



ducted, and quantitatively compared with various PN 
approximations. Quantitatively cross-checking waveform 
comparisons requires some care. How do the differences 
between PN and numerical results compare with intrin- 
sic numerical error estimate? How do different variants 
of PN predictions compare with the numerical results? 
Because each waveform comes with no meaningful abso- 
lute reference in time and phase, how can the freedom 
to offset these parameters be controlled for the compar- 
ison? Effective comparisons also require longer wave- 
forms than those produced in the earliest simulations. 
The first co mparison add r essing these issues came at the 
end of 2006 (jBaker et all 12007 dlh The study considered 
the case of equal-mass non-spinning black-hole mergers, 
comparing a 3.5 PN waveform with numerical results cov- 
ering the last seven orbits. The compared waveforms are 
shown in Fig. &)\ In the comparison the waveform phases 
agreed within 1 rad of phase drift, for a little over ten 
gravitational- wave cycles preceding the last orbit before 
merger, comparable to numerical error estimates. This 
result gave clear indication that PN waveforms could be 
accurate in the last orbits approaching merger, and that 
PN and NR, combined, could treat the complete wave- 
form signal. 

Not all equally valid variants (approximants) of the ap- 
proximate PN waveforms agree this well with numerical 
results. To understand this, it is worthwhile to consider 
PN results in a little more detail. The PN approxima- 
tion is formally understood as an expansion in powers 
of the speed of the merging black holes, v/c. A con- 
crete result of the theory will typically be a Taylor ex- 
pansion for a specific dependent variable, in terms of a 
chosen independent variable. The obvious choice, to ex- 
press the waveform itself (gravitational- wave strain) as a 
function of time, would give a poor result since the sinu- 



soidal shape of the waveform is difficult to approximate 
by a polynomial. PN waveform results are typically ex- 
pressed as separate expansions for the orbital phase, and 
polarization component amplitudes. These are given as 
expansions in the time to merger v/c cx t -1 / 8 or orbital 
frequency v/c cx fJ 1 / 3 . The orbital phase information 
may also be expresse d by an expansio n for Q, referred to 
as the chirp-rate; see iBlanchetl (2j006) for a review of the 
approach and a particular explicit waveform expansion. 
The PN information may also be encoded in an Hamil- 
tonian formulation description of the dynamics, which 
may be integrated numerically. Researchers also choose 
between Taylor series and other "resummed" expressions 
of the results based, f or instance, on Padc approximants 
(|Damour et all Il998l) . 

In the comparisons described above, the PN waveform 
phasing was found to agree with numerical results de- 
rived from a Taylor series expansion for chirprate 
expanded in powers of frequency. This is known as the 
TaylorT4 approximant in the language of the ground- 
based gravitational-wave community. Later studies with 
longer and more accurate numerical simulations provided 
more precise tests of phases and amplitudes for a variety 
of PN a pproximants (IBovle et all l2007al lHannam et all 
l2008d) . iBade et all (l2007al) found e ven closer agreement 



(than found in (|Baker et aUl2007dl) ) with the TaylorT4 
approximant, to within 0.05 rad over nearly 30 cycles be- 
fore the orbital frequency reaches MVl = 0.1 at roughly 
an orbit before merger. 

Figure [3D] shows the results of four PN approximants 
with terms kept to various PN orders. The excellent 
agreement of the TaylorT4 approximant with numerical 
results is not matched by other equally consistent PN 
waveform approximants in the late portions of the wave- 
forms. However, all four approximants agree to within 
about a radian if the comparison is cut off at waveform 



30 



frequency Muj = 0.1, about five orbits before merger 
[note that this is the frequency of the dominant (2, 2) 
mode, twice the orbital frequency MQ). 

While most attention has so far been given to the non- 
spinning cases, some results exist for other regions of pa- 
rameter space. For equal-mass configurations of black 
holes with spins aligned with the orbital angular momen- 
tum, the phase agreement of the TaylorT4 approximant is 
not generally best, and phase differences (over ten cycles 
preceding Muj = 0.1) are much larger than those seen 
in the nonspinning case, with d ifferences near 1 rador 
more for spins (a/M)i, 2 > 0.5 (|Hannam et all l2008ari . 
PN spin effects are not yet derived beyond 2.5 PN or- 
der, possibly limiting the accuracy of the PN waveforms 
in these comparisons. Some PN comparisons have also 
been examined for other configurations, including a pro- 
cessing system ( Campanelli et al\ . I2009T ) and for eccen- 
tric configura t ions o f equal-mass, non-spinning mergers 
([Hinder et aZ.l . l201fj ). 

In comparing Fourier-domain PN approximants with 
numerical results for a set of no n -spinn ing mergers with 
various mass ratios, IPan et al. I (120081) discovered that 
the abrupt truncation of the Fourier-domain waveforms, 
inducing Gibbs oscillations in the time-domain, can 
approximately emulate the physical quasinormal ring- 
down that terminates the numerical-simulation wave- 
forms. Even better agreement is possible if i], the sym- 
metric mass ratio (JTJ) , is allowed to extend beyond its 
physical range (to rj > 0.25) or if an extra pseudo-4PN 
order term is added. More recent studies have confirmed 
that these simple models are directly us eful in some 
gravi tational- wave observation applications ( Boyle et all 
l2009h . These results have motivated further development 
of phcnomenological full- waveform models. 



B. Analytic Full-Waveform Models 

For the early part of the waveforms various PN ap- 
proximants agree at 3.5 PN order, the best currently 
available, providing excellent approximations. At late 
times, the results of different approximants diverge, and 
the numerical-simulation results are the only way to accu- 
rately derive the predictions of general relativity. Can an 
analytic waveform description be developed that simulta- 
neously encodes the results of both complementary theo- 
retical approaches? The generally simple features of the 
numerically simulated merger waveforms raise hopes that 
fairly accurate approximate analytic waveform descrip- 
tions may be produced with simple dependence on the 
binary system parameters, providing an efficient means 
also of interpolating from a sparse sampling of the full 
parameter space for which numerical-simulation studies 
are completed. These waveform models would be compu- 
tationally efficient (compared to numerical simulations) 
and may have applications in a broad class of observa- 
tional data analysis applications. 

Ajith and collaborators have developed a Fourier- 



domain full- waveform model phenomen o logica l approach 
resembling the treatment of lPan et al. I (12008ft . but with 
greate r emphasis on matching the nu merical waveforms 
(lAiithl . 120081: lAiith et al. I. l2008ll2007ri . They be gin com- 
bining information from the PN and NR-based wave- 
forms by first stitching together time-series data for the 
dominant (£ = 2,m = 2) component waveforms, includ- 
ing a long PN precursor joined to a numerical merger 
waveform. Fourier transforms of waveforms are then 
fit to a parametrized model of the waveform similar to 
waveform families previously applied in ph enomenolog- 
ical t reatments of purely P N waveforms ( Arun et all 
120061 : [B uonanno et al. 1. 120031) . Their waveform model for 
mergers of non-spinning binaries has the form 
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where C is a constant related to the masses, and d is the 
distance to the source. The effective Fourier amplitude 
A e ff is modeled piecewise, 
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with distinct power-law segments before and after a 
merger frequency f m , and a Lorcntzian £(/, / r ,cr) de- 
cay beyond the transition frequency / r demarking the 
ringdown of the final black hole. The Fourier phase is 
represented by a single power-law expansion 



* off (/) = 2TTft Q + fo+J2*Pkf 



(fe-5)/3 



(38) 



For each waveform, the free amplitude and phase coef- 
ficients are determined by a fit to the stitched NR-PN 
hybrid waveforms. It was found that the mass-ratio 
parameter dependence could be modeled by simple fits 
to quadratic functions of the symmetric mass ratio (JlJ. 
These waveforms have been applied in gr avitational- wave 
data analysis stu dies (see Se c . IVII.AK Ajith and Bosd . 
120091 )). Recently lAiith et al\ (|2009t) have developed a 
generalization of this model for non-precessing system of 
spinning black holes. 

Stitching together PN and NR waveforms can be 
avoided, and closer contact with the basic physics main- 
tained, by requiring direct PN consistency in the full 
waveform model. Some PN approximants can resem- 
ble the numerical late-merger results. If such approx- 
imants can be tuned to agree with numerical results 
right up to merger by careful choice of adjustable pa- 
rameters which only affect PN consistency beyond known 
PN order, the result would simultaneously encode PN 
and numerical-relativity results. The effective-one-body 
(EOB) family of PN Hamiltonian models has promise 
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as a tu nable model for encod i ng; bo t h PN and NR 
resul t s dBuonanno and Damourl . 119991 120001: iDamouri . 
120081 : |d amour and Nagarl . l2010f ). In the time domain, 
techniques have also been developed for extending these 
waveforms into the ringdown of the final black hole. 
Analysis of the early numerical results in comparison 
with an untuned EOB model gave promising indications 
that the waveforms c o uld be closely approximated this 
way (jBuonanno et al I l2007al) . With tuning, it appears 
that this construction can provide an analytic but po- 
tentially very accurate, approximation to the complete 



coalescence waveform. 

In the EOB model, the binary motion is recast as 
the motion of a single effective body of mass /j = 
M1M2 / {Mi + M2) moving about a central potential, as is 
familiar from Newtonian mechanics. In the general rel- 
ativistic version of this framework, the effective body's 
motion follows a geodesic (to 2 PN order) around a mod- 
ified version of a Schwarzschild metric. The motion of 
the effective body is described by an effective Hamilto- 
nian, which, for systems of nonspinning black holes, may 
take the form 
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The expressions for A(r),D(r) and Q(r) are chosen so 
that the PN-expansion of the Hamiltonian is consistent 
with the results from PN t heory to known order (3 PN for 
nonspinning black holes) flBuonanno and Damourl . I1999L 
120001: I Damour et aD, I2000D . The nonconservative con- 
tribution to the motion, arising from the loss of angu- 
lar momentum to gravitational radiation, is encoded in 
an additional flux term entering as an external force in 
Hamilton's equations, which is also constrained to be 
consistent with PN theory (typically to 3.5PN order). 
Hamilton's equations are integrated to derive the effec- 
tive body's motion, and hence the black-hole trajecto- 
ries. Waveforms are constructed using PN extension of 
the quadrupole formula relating the motion of the black 
holes to the amplitude and phase of the multipolar radi- 
ation components. The last part of the radiation, arising 
after merger, is completed by continuously matching a 
sum of quasinormal ringdown modes to the waveforms 
truncated near the point of merger. 

While consistent with PN theory, the formalism can 
be adjusted to also match the wave forms derived from 
numerical-relativity simulations. The flux function, as 
well as A(r) and D(r), can be adjusted by the addition of 
higher-order PN terms to modify the strong-field dynam- 
ics to agree with NR without violating PN consistency. 
An early implementation of this approach, encoding the 
combined results PN and N R for nonspinning black - 
hole mergers was presented by iBuonanno et al. I (l2007bl) . 
Figure |3T] shows their result comparing numerical and 
adjusted EOB-model waveforms for a 4:1 mass-ratio 
merger. Subsequent work involving more accurate nu- 
merical simulations, and more careful tunin g of the EOB 
model has improved EOB model waveforms (iBovle et all 
l2008ri IDamour and Nagarl 120081: iDamour et all 120081) . 
Recent comparisons of improved EOB models with 
high- accuracy numerical results from lScheel et al. I (120091) 
yield differences comparable to numerical errors show- 
ing phase agreement within about 0.01 rad through 
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FIG. 31 Compa rison of EOB-model and numerical-relativity 
waveforms from IBuonanno et al\ (|2007bT l for the ft+ compo- 
nent of the gravitational-wave strain from the merger of bi- 
nary with mass ratio 4:1. The waveform is for an observer 
located an at inclination angle — n/3 from the axis of rota- 
tion. 



~ 30 gravitational-wave c ycles ( Buonanno et aZj . l2009"bt 
IDamour and Nagarl . [2009]) . 



Recently the 
to match cases 



EOB approach has been extended 

of non- precessing spinning black- hole 

mergers (|Pan et al ]. l2010h . There remain open questions 
in extending these waveform models to generic spin con- 
figurations, more extreme mass ratios, and full multipolar 
content. 




FIG. 32 Agreement of r efined EOB-model I — m = 2 waveform (|D amour and Nagarl . I2009T ) (red) with numerical result from 
Ref. (jScheel et all l2009h for equal -mass non-spinning merger (black). Slight differences near the merger time are difficult to 
perceive without color. 



C. Post-Newtonian Models for Numerical Initial Data. 

Before a numerical simulation can commence, a numer- 
ical relativist requires a model for the initial configura- 
tion of the two black holes and the geometric fields which 
represent them. As discussed in Sec. IIV.DI this involves 
not only solving the general relativistic constraint equa- 
tions, but also providing some ansatz for the free data. 
In some cases, some information from PN theory may 
be applied to produce an initial data model that more 
precisely represents the desired physical configuration. 

To understand how PN information can be applied, it 
will be useful to revisit the discussion of numerical initial 
data. Typical simulations begin with a modeling ansatz, 
the Brandt-Brugmann model, for instance, which pro- 
duces initial field data from a given set of particle-like 
parameters related to the black holes' masses, spins, ini- 
tial positions, and momenta. Usually, the initial data 
models make no attempt to represent the gravitational 
radiation fields which would have been previously gener- 
ated by the motion of the black holes, and do not uti- 
lize PN information in this reduction of the field degrees 
of freedom to particle parameters. Information from PN 
theory is frequently applied in choosing the specific parti- 
cle parameters which correspond to the sought-after sim- 
ulation, particularly for circular inspiral configurations. 
Work has also begun toward making richer use of PN 
information for improving the field ansatz. 

Most of the simulations discussed have been designed 
to represent black holes in circularized orbits. Before 
long-lasting numerical-simulation results were available, 
comparative studies of numerical initial data sets with 
PN-derived informat ion provided a gauge of the re- 
sults [see for instance iBaker et al. I (|2002bh :T Caudill et al\ 
(120061) ; ICook and Pfeifferl (|2004l) ]. Without evolutions, 



can lead to eccentricity in the simulatio n which im- 
pacts t he simulated radiat ion waveforms ( Boyle et all 
l2008bl: iPfeiffer et all 120071) . The residual eccentricity 
can be reduced via a straightforward iterative procedure 
(|Pfeiffer et aU l2007lh but this requires repeated simula- 
tions lasting several orbits; expensive in both time and 
computational resources. To minimize this eccentricity 
without resorting to iteration, several groups use trajec- 
tory information from PN theory in setting up the pa- 
rameters for the numerical initial data. This approach 
enables simulations with eccentricities e < 0.002 for 
equal-mass n on-spinning merger s with initial separations 
di„it > 10M (jHusa et aq|2008b| ). The technique is help- 



these studies focused on theoretical properties of the 
black-hole configuration space such as the ISCO (see 
Sec.lILATJ). 

Still subtle imbalances in the initial data, such as 
excess angular momentum for the chosen separation, 



ful in simulations with unequal masses and nonvanishing 
spins a s well, though the res idual eccentricity is generally 
larger (jWalther et all 120091) . 

Work is underway on techniques allowing the richer ap- 
plication of PN information from numerical initial data 
modeling. The common assumption of a conformally flat 
spatial metric in numerical initial data, disagree s with 
PN results at the 2 PN order (jDamour et aUl2000t) . mak- 
ing this a likely leading source of modeling error which 
produces spurious initial transients in the simulations. 
Though not yet as well-developed as the widely applied 
models, an alternative approach applies metric informa- 
tion fro m PN theory for a non - trivial initial con f ormal 
metric ( Jaranowski and Schaferl . Il998l ; lOhta et all . 11974 
ISchaferl . 1985 ). With these techniques it is also possible 
to encode in the initial data information about the prior 
radiation generated by the system before the "initial" 
time of the simulation dKellv et all l2007t iTichv et all 
I2003D . 



D. Post-Newtonian Theory for Interpretation of Numerical 
Results. 

We have noted several specific areas where research 
in numerical simulations makes contact with PN theory. 
These alone do not provide a full picture of the interplay 
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between the two theoretical approaches. Most broadly, 
PN theory provides a foundation for interpreting numer- 
ical results. 

Numerical relativists draw widely from PN-based back- 
ground knowledge of the black-hole binary. Many impor- 
tant phenomena, such as the ISCO, spin-orbit coupling, 
spin precession, and 'gravitational rocket' kicks were first 
studied in the PN approximation, providing a foundation 
for subsequent numerical studies. Small examples can be 
noted throughout this review. In particular we point to 
the value of PN results in interpreting mergers of spin- 
ning black holes (Sec. IV.D|) and kicks (Sec. IV. C[) . Ana- 
lytic formulas for approximately expressing the final state 
of the black hole, mass, spin, and momentum have also 
drawn heavily from insights based on the PN treatment 
(see Sec JVIILA.il ) 



VII. APPLICATIONS TO GRAVITATIONAL WAVE DATA 
ANALYSIS 

Several detectors are active or in their planning stages 
to detect the gravitational-wave signals from astro- 
physical processes. Prominent among these are the 
ground-based i nterfero metric detectors - LIGO (U.S.A.) 
(j Abbott et all , l2009al ). GEO (Germany), Virgo (Italy), 
TAMA (Japan), and AIGO (Australia) - sensitive to 
frequencies in the range 10 1 — 10 3 Hz, as well as next- 
generation instrum ents such as the Einstein Telescope 
( Freise et all l2009h . Also planned to launch at the end 
of the next decade is the Laser Interferometer Space An- 
tenna (LISA), with complementary frequency sensitivi- 
ties between 10 -4 — 10 _1 Hz. All of these instruments 
are subject to a variety of noise sources; in the ground- 
based detectors, these sources will completely overwhelm 
the signal, unless filtered intelligently. Figure [33J shows 
the design sensitivity curve resulting from these disparate 
noise sources, for the LIGO and Virgo detectors. 

Burst analysis can extract many signals, but op- 
timal results in gravitational-wave data analysis re- 
quire matched filtering of t he noisy input signal 



(iFlanagan and Hughes! . Il998l IWainstein and Zubakovl 
Il962h . Fundamental to this approach is some measure 
of overlap between a physical signal h±(t) and filtering 
template waveform h 2 (t) . For this it is convenient to use 
a freque ncy-space inner product (-|-) between signals, de- 
fined as ( Cutler and Flanagan! . Il994l) 
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where hi(f) and h 2 {f) are the Fourier transforms of the 
signals, and S n (f) is the (one-sided) noise power spectral 
density of the detector we are interested in. This can be 
norm alized to produce the match (or overlap) (jOwenl . 
I1996D between two waveforms. 
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FIG. 33 Design sensitivity curves for the LIGO & Virgo de- 
tectors (dashed lines), as well as the approximate curves used 
for the NINJA project (solid lines). LIGO sensitivity is ef- 
fectively zero below ~ 30Hz, while Virgo does a little better 
(note that achieved sensitiv ities may not perfec tly mirror the 
design curves). Taken from lAvlott et all (|2009h . Reproduced 
by permission of the Institute of Physics. 



To use (|4"0")) requires the availability of a set of wave- 
form templates - simple, few-parameter analytic model 
waveforms to filter against. Even for a relatively sim- 
ple system such as a black-hole binary, the combination 
of possible source and detector configurations leads to 
a 17-dimensional parameter space 5 . Some of these pa- 
rameters are intrinsic to the source: total binary mass 
M = Mi + M 2 , symmetric mass ratio r\ ([1]), spins Si, 2, 
time to coalescence t c , eccentricity e, and eccentric phase 
<f) e . The remaining parameters have to do with the rela- 
tion between source and detector: (luminosity) distance 
to source Dl, inclination t, orbital phase <f>, waveform 
polarization ip, and position of the binary on the detec- 
tor's sky {0,$}. Restricting consideration to binaries 
that have circ ularized by th e time they enter the detec- 
tor's window ([Peters! . Il96l allows us to neglect the ec- 
centricity parameters {e, 4> e }. Crucially also all observ- 
ables have a simple dependence on the total (rcdshiftcd) 
mass, (l+z)M. Similarly, the observed waveforms have a 
trivial dependence on the time to coalescence t c . As this 
means only one theoretical waveform has to be generated 
to cover all astrophysical masses and coalescence times, 



5 It is easy to see that there are a total of 17 dimensions. Gen- 
erally, the instantaneous state of each black hole has ten de- 
grees of freedom for its mass, position, momentum and spin. Of 
the 20 degrees of freedom for two black holes, the three related 
to the center-of-mass momentum are unmeasurable by distant 
gravitational- wave observations. The peculiar motion has neg- 
ligible effect, and the proper motion results in a Doppler shift 
indistinguishable from a change in the total mass. Also note 
that changes in eight more degrees of freedom, corresponding to 
the center-of-mass position in spacetime, spatial orientation and 
total mass have trivial effects which can be treated analytically, 
leaving, in full generality, a space of nine parameters that must 
be covered by simulations. 
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(l+z)M and t c are sometimes treated as extrinsic param- 
eters instead. Nevertheless, numerical relativists are left 
with a seven-dimensional parameter space to cover (77, 
Si, S2), with no obvious shortcuts to lighten the work- 
load further. Moreover, numerical methods will always 
struggle with finite accuracy. 

Equipped with a parametrized set of template wave- 
forms h m (t;X), we can test an incoming detector data 
stream s(t) = h e (t) + n(t), consisting of a (pos- 
sible) "exact" gravitational wave h e (t) and detector 
noise n(t) by calculating the signal-to-noise ratio (SNR) 
(jCutler and Flanaeanl [l994h : 



yj (h m (\i)\h m (\i)) 



(41) 



In the case where the signal h e (t) in the data stream 
corresponds perfectly (up to overall scaling) with the 
model waveform h m (t; Xi) for some value of the model 
parameters Xi = A^, we obtain the optimal SNR: p opt = 
\J{h e \h e ). We refer to this simply as the "SNR" p for 
the remainder of this section. 

Expected SNRs depend strongly on the detector in 
question. Generally speaking, the ground-based detec- 
tors in operation (LIGO, GEO, VIRGO, etc) are ex- 
pected to observe with only modest optimal SNRs, while 
LISA is ex pected to observe MBH mergers with SNRs of 
hundreds (jDanzmann et all , Il998l ) . 

Of course, real model waveforms will be inaccurate, 
due to incomplete knowledge of the underlying physics, or 
perhaps a desire for simp licity. In this case, the achieved 
SNR cannot be optimal. lApostolatosI ( 19951 ) denned the 
fitting factor (FF) as the reduction in signal from such 
imperfect templates: 
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A fitting factor of 1 means that the exact signal lies 
somewhere in the model space. However, for a particular 
detector, it is impossible to generate templates arbitrarily 
close together in any one parameter dimension. We can 
develop a sense of how closely spaced a set of templates 
must be so that any physical signal will b e detected by 
at least one of them with high likelihood. lOwenl ( 19961 ) 



developed a geometric picture of a template-space metric 
defining the distance between neighboring closely-spaced 
templates. With these tools, he was able to link the 
number of templates Af to the desired minimal match 
(MM) - the worst value of the match between the signal 
and any template. The c onnection between thes e con- 
cepts is nicel y laid out by Lindblom et all (2008). Fig- 
ure [SI from ILindblom et al\ ~i 20081 ). demonstrates how 
£mm = 1 — MM and £ff = 1 — FF combine to produce 
an effective mismatch £eff- 

The percentage of detected signals scales as the cube 
of the fitting factor. For this reason, detector ana- 
lysts require rather stringent fitting factors from their 




FIG. 34 Relationship between exact waveform h e , model 
waveform with the same physical parameters h m , actual "best 
fit" model waveform hf h, and template bank waveforms ht, 
hy . [Figure taken from ILindblom et al] (|2008l )] 



template banks, e.g., 97% (implying detection of 90% 
of all signals pre s ent). For initial LIGO, for instance, 
iBaumgarte et al. I (120081) estimate that ~ 100 zero-spin 
templates (produced from ~ 10 different numerical sim- 
ulations) would suffice to detect all nonspinning merg- 
ing binaries. Obviously we expect that the addition of 
spins of significant magnitude and arbitrary direction will 
greatly increase this number. 

Detection is only a first step, although an important 
one. When signal strengths are high, as is expected with 
LISA for massive binaries, we might expect to be able 
to read off some of the system's intrinsic and extrinsic 
parameters. Achieving this requires templates that ac- 
curately cover the parameter space of the binary, with 
a one-to-one relationship with the parameters of the bi- 
nary - that is, that the best-fit model waveform hfh and 
closest-parameters model waveform h m in Fig. [33] coin- 
cide closely (within statistical errors). It also requires 
that we have an understanding of the probabilistic cor- 
relations between parameters. 

For high-p signals, the simplest way to estimate param- 
eter uncertainties is using the Fisher inf ormation matrix 
( Cutler and Flanaeanl . Il99l iFinnl . Il99l . which uses sim- 
ple derivatives of the template waveform shape with re- 
spect to the physical parameters: 
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In the high-SNR limit, the parameters are assumed to 
have a multivariate Gaussian distribution, and the Fisher 
matrix is the inverse of the covariance matrix between pa- 
rameters. More sophisticated approaches that do not as- 
sume high SNR include Markov-chain Monte Carlo meth- 
ods dGilks et aZ.l . ll996l) . 

Different requirements can be made on the quality of 
the de veloped templates, as described by iDamour et al\ 
( 19981 ). The first is that they be effectual. Roughly 
speaking, this means that the templates must be capable 
in the bulk of detecting gravitational- wave signals. In 
the language of matched filtering above, we demand that 
the fitting factor FF be extremely close to unity for any 
signal. 

More stringently, we might demand that the templates 
be faithful: that is, that the template waveform cor- 
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responds to underlying source parameters that corre- 
spond closely to the parameter s of the detected wave - 
form source. In the language of iLindblom et all ( 2008j ). 
the model waveforms h m and hfh coincide closely (see 
Fig. I34p . This requirement will be crucial for identifying 
the astrophysical source of the radiation. 

Traditionally, templates used in the search for 
comparable-mass black-hole binaries have depended 
wholly on PN theory. As detection strategies arc most 
sensitive to phase discrepancies between template and 
observed signal (see below), the determination of wave- 
form phase to high accuracy has been crucial. Modern 
PN templates include phase terms up to 3.5PN order for 
nonspinning binaries, with 2.5 PN order spin corrections. 
Many templates use this high-accuracy phase with the 
leading-order quadrupole amplitude to produce restricted 
templates; however, amplitude corrections are k nown to 
2.5 P N order (higher for certain \ow-£ modes) (|Kidderl . 
I2008h . 

While such PN-driven templates offer excellent phase 
accuracy during the long, slow inspiral of the binary, their 
usefulness becomes questionable as we approach the last 
pre-merger orbits of the binary. We may ask, then, what 
does our new numerical insight into the final moments of 
merger bring us? 



A. The Direct Impact of Merger Waveforms in Data 
Analysis 



Once full numerical waveforms became available, sev- 
eral groups used the results to test older predictions of 
the effect of the merger segment on observability. For 
the equal-mass nonspinnin g case, the numerical merger 
gave results consistent with lFlanagan and Hughesl (|l998l ) 
for initial LIGO, though with a smaller m erger SNR 
dBaker et all l2007d: iBuonanno et all l2007a|) . The boost 
in SNR from merger is significantly greater f or Advanced 
LIGO and LISA, as shown in Fig. [35j from lBaker et al] 
(l2007dh 

The relevance of the final plunge and merger to the 
overall SNR of the binary depends on the total mass of 
the binary and where it falls in the window of th e detec- 
tor. Figure [31)1 adapted from lBaker et al. I (l2007d) . shows 
this for LISA, plotting the "characteristic signal strain" 
/i c har(/) = 2/| h(f) | (for the dominant quadrupole radia- 
tion). Lower-mass binaries have higher- frequency wave- 
forms at all dynamical stages: for an equal-mass binary 
with M < 1O 4 M0, the late- merger and ringdown signal 
will fall outside LISA's sensitivity band. For such low 
masses, inspiral-only waveforms should be adequate for 
data-analysis purposes. 

The LIGO and Virgo detectors are sensitive to fre- 
quencies above ~ 30Hz (see Fig. [33]). The merger 
of binary systems is marked by a chirp gravitational- 
wave signal, whose monotonically increasing frequency 
saturates at the dominant quasinormal mode, which de- 
pends on the post-merger hole's mass and spin. For 
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FIG. 35 The importance of including NR merger waveforms 
in data analysis can be seen h ere for Advanced LI GO (top) 
and LISA (bottom), taken from lBaker et al\ (|2007d ). In both 
panels, the dashed curve is the achieved SNR (for a given 
rcdshifted source mass) when only the early inspiral signal (up 
to ~ 1000M before merger) is used; the dotted curve uses the 
signal between this time and when final plunge commences 
(1000M to 50M before merger; the thin solid curve uses just 
the merger waveform (starting 50M before merger). The thick 
solid curve is the combined result of using the entire waveform. 
The SNR, and hence distance reach, of the detector is clearly 
greatly enhanced by the final merger portion. 



a nonspinning binary of total mass M, this /qnm ~ 
1.75 x 10 4 (Af/M Q ) _1 Hz. This means that a merging bi- 
nary of mass greater than ~ 600M Q will never be seen by 
LIGO. In contrast, most post-Newtonian templates stop 
at /iSCOj the frequency at the innermost stable circular 
orbit (well defined only for test particles; see Sec. M.ATTj) . 
This is 



/isco 



1 



6V6M 



1.36 x 10 4 (M/M o )- 1 Hz. (44) 



For an incomplete PN waveform to be useful, we want 
the missing merger-ringdown section (i.e. the part of the 
signal with / > /isco) to contribute as little as possi- 
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FIG. 36 The importance of binary mass for the LISA ob- 
servati on window, demons trated by the characteristic ampli- 
tudes (jBaker et al I l2007ch of three different sources, relative 
to the rms noise amplitude of the LISA detector. On each 
hchar curve, we mark times before the peak amplitude (cir- 
cle) - one hour (square), one day (diamond) and one month 
(triangle). 



ble to the total SNR. For instance, if we take the high 
end of LIGO's sensitivity window to be ~ 800Hz, this is 
/rsco for a binary of total mass M ps 17Mq. Then, for 
binaries with M > 17 M & , f > /isco will fall within the 
visible window; that is, LIGO will see the final stages of 
the binary merger, where no acceptable PN wavefo rm is 
available. Recent work by iBuonanno et al. | (l2009al) sup- 
plied a more precise answer: inspiral-only PN templates 
can be trusted for all LIGO configurations for systems 
with M < 12M Q . Above this critical mass, there is a 
gap in reliable information, which must be filled by nu- 
merical results. 

Figure [371 adapted from |Pan et al. I (I2008D . shows the 
importance of the last stages of merger in initial LIGO 
as a function of mass range. The thin blue curves are hy- 
brid NR/EOB waveforms, while the thick red waveforms 
are generated from these by "whitening": that is, they 
have been Fourier-transformed to the frequency domain, 
then rescaled by 1/ \J S n (f), and finally re-transformed 
to the time domain. The whitened amplitude in a seg- 
ment is proportional to the contribution of that seg- 
ment to the total SNR; each marked segment accounts 
for 10% of the total. As the total system mass M in- 
creases, the whitened amplitude becomes more bulked to- 
ward the merger time. For the largest total mass shown, 
M = IOOM0, more than 90% of the signal power comes 
from the last cycle + merger + ringdown. 

In Se c. IV.B.31 we reference d the recent "Samurai" 
project ( Hannam et all l2009al ). which establishes the 



consistency of all of the "long" (pre-merger duration 
•> 1000M) equal-mass nonspinning waveforms within 
their stated numerical accuracy. As well as direct com- 
parisons of phase and amplitude errors over time, they 
conducted mismatch tests in the regime of the LIGO and 
Virgo ground-based detectors. They demonstrated that 
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FIG. 37 The power distribution in 10% segments of the wave- 
form before merger, for a range of LIGO-appropriate equal- 
mass binaries. Note that as the total mass M increases, the 
final merger and ringdown acc ount for a greate r fraction of 
the total power. Adapted from lPan etoH (|2008T l. 



quadrupole waveforms (£ = 2,?n = ±2) from the five 
numerical codes have tiny mismatches (~ 10~ 3 ) for bi- 
nary masses above 6OM with "enhanced" LIGO, and 
masses above 18OM with Advanced LIGO, Virgo and 
Advanced Virgo. All waveforms would be indistinguish- 
able for SNRs below 14, entirely reasonable for the cur- 
rent generation of ground-based detectors; this is shown 
for the Advanced LIGO detector in Fig. [3HJ 

Studies using incomplete (non-merger) PN waveforms 
have shown that introducing fuller harmonics of spin- 
ning and precessing binary systems greatly enhances pa- 
rameter estimation, removing some parameter correla- 
tions and hence reducing certain errors by several or- 
ders of magnitude (tho ugh other parameters are sig- 
nificantly less improved) (lLang and Hughes . 20061 20081 



120091 : ISintes and Vecchiol l200Ct iTrias and Sintesl . 120081) 
It is expected that the introduction of full numerical 
waveforms for the merger portion will lead to further 
improvements, at least in regions of the detector's win- 
do w where the binary m erger is visible. A recent study 
bv lAiith and Bosel ( 2009T ) used phenomenological merger- 
inspiral-ringdown templates (discussed in Sec. IVII.Bp to 
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FIG. 38 Maximum signal-to-noise ratio (SNR) p below which 
different "long" numerical waveforms are indistinguishable 
for the Advanced LIGO detector. Even the earliest, lowest- 
accuracy published long waveforms are indistinguishable for 
SNR levels belo w 14, as demonstrated in this figure from the 
samurai paper (|Hannam et al.l,l2009ar i. 



assess these improvements for Advanced LIGO and Ad- 
vanced Virgo. They showed that parameter-estimate ac- 
curacy is improved substantially compared to inspiral- 
only waveforms for system masses > 20M Q ; in particular, 
the average sky-position error for an M = lOOM© equal- 
mass, nonspinning binary at a luminosity distance of 1 
Gpc drops to about one-tenth of a square degree. Since 
they only uses the dominant mode for this investigation, 
there will be systematic errors as well. 

Numerical verification of the extent of parameter es- 
timation i mprovements for LISA has led to differing 
estimates (iBabak et all 120081 : IMcWilliams et all l20irj 
iThorpe et aZj . l2009h . In particular. iBabak et all (|2008h 



find that 50% of injected signals could be located on the 
detector sky within an error box of 3 arc min, while t c 
could be measured to less than a second and Dl within 
1.5% (ignoring issues of weak lensing). They attributed 
these impressive results to a combination of very high 
SNR (between 900 and 9000) and the use of higher mul- 
tipole s of the radiation. In contrast, IMcWilliams et al\ 
( 20101 ) found smaller (by about one order of magnitude) 
improvements over inspiral-only waveforms. This may in 
part be attributable to the use of different mass ratios 
and different-length waveforms by the two groups, but 
has not been fully resolved. 



B. Developing Analytic Inspiral-Merger-Ringdown 
Gravitational Waveform Templates 

One of the eventual aims of numerical-relativity sim- 
ulations of binary mergers is the production of a set of 
gravitational-wave "templates" that will cover the many- 
dimensional space of astrophysical parameters. Numeri- 
cal merger simulations are still orders of magnitude too 
slow to run on the fly to filter incoming detector data 



streams; thus there is a need to develop analytic wave- 
form expressions that encode the numerics. 

It has been demonstrated i n the physical regime o f 
nonspinning binary mergers (|Buonanno et all l2009ah 
that inspiral-only PN-based templates become inconsis- 
tent with each other (and with full waveforms) at to- 
tal binary masses M > 12M Q for initial and advanced 
LIGO configurations; above this mass, full templates in- 
cluding the numerical-simulation-based understanding of 
merger / ringdown become necessary. 

One such set of tem plates incorporati ng numerical- 
relativity data is due to iPan et all (l2008f ) . These were 
extensions of the Stationary-Phase Approximation (SPA) 
templates used in LIGO/ Virgo data analysis, with input 
from equa l-mass numerical wavefo rms supplied by Frans 
Pretorius ( Buonanno et all\2007e$) , as well as equal- and 
unequal-mass numerical waveforms from the Goddard 
group. These templates are effectual - they achieved 
fitting factors of > 0.96 - but they did this by extending 
source parameters into unphysical space. 

To date a handful of useful faithful full-waveform 
template banks has been produced, aimed at covering 
all mass ratios for nonspinn ing binaries . The "phe- 
nomenological" templates of lAiith et all (|2008l . l2007h 
are simple three-segment curves in frequency space, 
with matching paramet ers tuned by numerical data 
generate d by the BAM dBriigmann et all l2008b|) and 
CCATIE dPollnev et all l2007h codes. The te mplates of 
iBovle et all (|2008bf ); iBuonanno et ail (|2007bl ). coded in 
the LSC Algorithm Library as "EOBNR", are an ex- 
tension of cffcctivc-one-body waveforms, with a single 
"pscudo-4PN" parameter tuned by numerical data (ini- 
tially from the Goddard Hahndol code, and later using 
waveforms from the Caltech-Cornell group. Both the 
phenomenological and EOBNR template banks are dis- 
cussed in Sec. IVI.BI 

Both template banks are faithful in the sense de- 
scribed above, in the restricted (M, rf) parameter space, 
and both have been used in data-analysis in jection tests 
(jAvlott et aU l2009t ISajitamaria et all . l2009() . Currently, 
EOBNR templates are being used both for injection and 
filtering in the high-mass region of LIGO test analy- 
sis, and also for injection for the ringdown band. They 
were also both used in the matched-filter analyses of the 
NINJA project (see Sec. IVII.CI) . performing as well as 
the inspiral-only templates in detection, and considerably 
better in the li mited param e ter es timation attempted. 
Figure [331 from lAvlott et all ( 20091) . shows the accuracy 
with which the total mass M and time of merger were 
extracted from all detected injections in NINJA, when 
the EOBNR templates were used. Results for the phe- 
nomenological templates were similar for M, but per- 
formed slightly less well for time of merger, as it is not 
an explicit part of this model (given the relatively small 
number of samples in the NINJA tests, such minor dif- 
ferences may not be significant). 

As significant spin is expected in astrophysical black 
holes at all scales, it is important to develop templates 
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FIG. 39 Accuracy in parameter estimation for EOBNR 
templates in the NINJA project. The top panel shows 
the fractional error in estimated total mass, (Mi n j cctc( j — 
A^detected)/Afinjected, for all detected signals, w hile the lower 
shows the error in merger time. Adapted from lAvlott et al\ 
(|2009h . Reproduced by permission of the Institute of Physics. 



that encode spin also. IVaishnav et all ( 2007t ) demon- 
strated that matched filtering with nonspinning merger 
waveforms does a good job in detecting waveforms from 
binaries with significant spins. As their spinning wave- 
forms come from aligned-spin mergers where the total 
spin is zero (one hole's spin is aligned with the orbital an- 
gular momentum, while the other's is anti-aligned) , their 
configurations will lack some prominent spin-related ef- 
fects, including spin-orbit effects (pulling in or pushing 
out the ISCO of the orbit) and precession of the orbital 
plane. Their conclusions, then, are likely to be opti- 
mistic, and more generic configurations can be expected 
to do much worse with nonspinning templates. The phe- 
nomenological templates of Ajith et al. have recently 
been extended to include non-precessing spinning sys- 
tems (|Aiithl . 120091) . These have been shown to be both 
effectual and faithful in searches over non-precessing bi- 
nary hybrid waveforms with LIGO, where earlier zero- 
spin templates failed badly. For maximum simplicity, 
these new phenomenological template families are ulti- 
mately parametrized by three numbers: the total mass 
M, the symmetric mass ratio rj |T]), and a single effective 



{a/MU = £-±^(a/M)i + ^-^-(a/M) 2 , (45) 



where S = (M x - M 2 )/M = y/l - 4?? 2 



C. Using Numerical Waveforms in Data Analysis 
Applications 

An obvious application of numerical waveforms for 
data analysis is in directly testing analysis algorithms. As 
current templates are overwhelmingly based on PN infor- 
mation - either extended heuristically to include merger 
and ringdown, or terminating before merger - so current 
detection and parameter estimation algorithms are also 
based on incomplete information. 

Recently, a multi-group collaboration (the NINJA 
project) has performed the first direct insertion of ex- 
plicit numerical waveforms into a simulated LIGO datas- 
tream, and investigated detectio n efficiency and system- 
atic errors ( Avlott et all , 120091 ). Essentially all active 
numerical-relativity groups contributed waveforms; these 
were scaled randomly to represent different-mass sources, 
resulting in 126 different injections. Nine data-analysis 
groups analyzed the post-injection data stream, with a 
variety of methods. These fell into three main groups: 
matched filtering, burst analysis, and Bayesian tech- 
niques. 

For the matched filtering, some groups used stan- 
dard "TaylorF2" inspiral-only templates generated by 
PN theory, some used ringdown-only templates, and 
some used full inspiral-mcrgcr-ringdown templates, such 
as the EOBNR and phenomenological waveforms dis- 
cussed below and in Sec. IVI.B1 Both the TaylorF2 and 
the full template waveforms yielded quite high detection 
rates - around 80 of the 126 injected signals in triple coin- 
cidence between the LIGO detectors - while the ringdown 
templates performed more poorly, detecting only 45 sig- 
nals in triple coincidence. However for parameter estima- 
tion of the total mass M, most of the TaylorF2 estimates 
have an error of 40% or more, significantly more than 
for the EOBNR and phenomenological estimates. Fig- 
ure [35] shows errors in detected mass (upper panel) and 
coalescence time (lower panel) when using the EOBNR 
templates. While mass errors are still quite large, be- 
tween ~ —30% and +50% for the bulk of detected sig- 
nals, they are significantly less biased than t hose from the 
Taylo rF2 templates (compare with Fig. 8 of lAvlott et all 
(2009)). Estimation of the coalescence time t c was better 
for the EOBNR and phenomenological template banks 
as well, as might be expected (since TaylorF2 templates 
must necessarily cut off before coalescence). However, 
the EOBNR and phenomenological templates typically 
covered a larger mass range as well (up to 200 Mq com- 
pared to 9OM for TaylorF2), which will bias the results 
when the injected signals had a large mass. 
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Burst analysis is a more abstract approach to the 
problem, using partial information about signals, such 
as dominant frequency and approximate duration. In 
NINJA, two burst techniques were used - the "Q 
pipeline" and the Hilbert-Huang transform (HHT). 
The Q pipeline fits signals to a sine-Gaussian model 
parametrized by a central frequency and a Q factor 
(~ /r, the number of oscillations under the Gaussian). 
At the single-detector level, the detection rate is com- 
parable to that of the EOBNR and phenomenological 
matched-filter searches. Parameter estimation is limited 
to estimates of the peak power time (close to t c ) and peak 
frequency; for these parameters, performance was compa- 
rable to the ringdown matched-filter search. The Hilbert- 
Huang transform decomposes input data into "intrinsic 
mode functions" , each characterized by a single frequency 
scale. Applied to the signal, the HHT gen erates a high- 
resolution time-freq u ency map of the data ( Camp et all 
120071 ; IStroeer et al\ . 120091) . The HHT applied to the 
NINJA data found ~ 80 coincident events, competitive 
with the best matched-filter results; however, no param- 
eter estimation is currently possible with this method. 



Bayesian analysis attempts to reconstruct the poste- 
rior probability density function of a parametrized signal 
based on the observed signal. Two variant Bayesian ap- 
proaches were taken in NINJA analysis, involving dif- 
ferent Monte Carlo methods - Markov Chain Monte 
Carlo (MCMC) and Nested Sampling. In both cases, a 
parametrized waveform model is needed. For the MCMC 
study, a restricted PN waveform was used, only 1.5 PN 
in phase, but including leading spin effects. Detection of 
injected signals performed well, but the high mass of the 
injections meant that the bulk of the SNR occurred in 
the merger-ringdown phase; as a result the inspiral-only 
model masses merger time was significantly biased. The 
Nested Sampling effort instead uses the TaylorF2 tem- 
plates (restricted amplitude, 2.0 PN in phase) as well 
as phenomenological templates. In detection, the phe- 
nomenological templates significantly out-perform the 
TaylorF2 model; in parameter estimation, it achieves re- 
sults consistent with those of the matched-filter analysis. 



Overall, current results are encouraging, with high de- 
tection rates of injected waveforms, and good, if limited, 
accuracy in estimated parameters when attempted with 
faithful template banks. However, the project was ham- 
pered by the limited length of the contributed waveforms 
and unrealistically stationary LIGO and Virgo noise pro- 
files. Follow-on work in this area will demand longer nu- 
merical waveforms that span the LIGO sensitivity win- 
dow, or more likely, a systematic blending of long PN in- 
spiral waveforms to late-inspiral/mcrger numerical wave- 
forms, but also may be extended to other gravitational- 
wave sources. 



VIII. IMPACT ON ASTROPHYSICS 

The recent successes in modeling of binary black-hole 
mergers have captured the attention of astronomers and 
astrophysicists. While black-hole binaries have been 
studied both observationally and theoretically for many 
years, most efforts have been primarily within the frame- 
work of Newtonian gravity. The new developments in 
numerical relativity now supply the missing piece: the 
effects of the final merger in the strong-field general rel- 
ativistic regime. In this section we highlight three key 
areas of current interest: recoiling black holes, the spin 
of the merged black hole, and electromagnetic signatures 
from the final merger. 



A. Recoiling Black Holes 

One of the more dramatic implications of gravitational- 
wave computations for astrophysical observations is the 
possibility that gravitational radiation-induced recoil will 
eject black holes from their host galaxies. The largest 
recoil velocities predicted by numerical relativity exceed 
the escape velocities of many galaxies (~ 1000 km s" 1 ). 
Calculating the probability of this kind of event requires a 
detailed understanding of the recoil on mass ratios, spin 
magnitudes, and spin orientations, together with some 
expectations about the distributions of these parameters. 



1. Predicting the Recoil 

It would be useful to have in hand a simple, analytic 
formula expressing the relevant dependencies of the re- 
coil velocity. The highly nonlinear, strong-field interac- 
tions of merger dynamics determine the final recoil ve- 
locity of the merged remnant. Numerical simulations 
are required to accurately compute such results. How- 
ever, in attempting to construct an ansatz for a phe- 
nomenological formula, one might initially assume a form 
consistent with PN predictions. This will suggest the 
functional dependence on mass and spin, and their lead - 
mg powers. For example, PN analysis (jFitchettl . Il983f) 
suggests that the recoil due to unequal masses may be 
proportional to v = Arj 2 y/1 — 477(1 + Br/) (the "Fitch- 
ctt formula"), where r\ is the symmetric mass ratio (fT|. 
Meanwhile, to leading PN orders, the spins should con- 
tribute through components of A = S2/M2 — Si /Mi and 
S = Si + S? (ICampanelli et all l2007bl iKidderi . 119951 
iRacine et qfl . l2009h . 

Such an ansatz can be further supported and con- 
strained by symmetry arguments. An economical ap- 
proach is to consider the recoil velocity as a Taylor expan- 
sion in all six spin components (three for each black hole) , 
where the coefficients are functions of mass and initial 
separation. Rotation, parity, and exchange symmetries 
of the binary syst em then impose relationships between 
these coefficients ( Boyle and Kesdenl 120081 : iBovle et all 
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l2008al ). By such arguments, for example, the compo- 
nent of the recoil velocity parallel to the orbital angular 
momentum can be shown to be proportional, to leading 
order in spin, to the dot product of the coordinate separa- 
tion vector between th e black holes a nd some linear com- 
binati on of the spins ([Baker et al I l2008at IBovle et all . 
l2008aTl . 

Unknown coefficients in the ansatz can then be 
fit to numerical results, and correction terms can be 
added as needed for numerical agreement. With these 
considerations, a tentative formula has been taking 
shape through the combined efforts of th e numerical 
relativity community, of the following form dBaker et al 



2007 



. I2008al 

2007d: iGonzalez et all. l2007bl: lLousto et all. I201C ; 



iBriigmann et al\, 2008a ; Campanelli et al. 



Lousto and Zlochowerl l2009t Ivan Meter et al. . 2010a ): 



K-ocoii = v m e\ + vj_(cos£ei + sin£e 2 ) + U|| e 2 , (46) 

(47) 



v ± = tl- r 



T ( || || 



(1 + 9) 
K 2 y 2 + K 3 if 
9 + 1 



(48) 



x [qa^ cos(0i - $i) - cos(0 2 - $2)] 
K s (g-l)r) 2 
(q + I) 3 

x [q 2 a^ cos(0i - $1) + a? cos(0 2 - $ 2 )](49) 



where v m is the contribution due to mass asymmetry, v± 
is the contribution due to spin that yields a kick per- 
pendicular to the orbital angular momentum, i>|| is the 
contribution due to spin that yields a kick parallel to the 
orbital angular momentum, oq is the projection of the di- 
mensionless spin vector d?i = Si /Mf of black hole i along 
the orbital angular momentum, af- is the magnitude of 
its projection into the orbital plane, 4>i refers to the 
angle made by af- with respect to some reference angle 
in the orbital plane, and $1 and $2 are constants for a 
given mass ratio and initial separation. Here the spins 
are considered to have been measured at some point be- 
fore merger or ideally, arbitrarily close to merger. $1 and 
$2 encode the amount of precession of each spin before 
merger. 

One can then evaluate the formula over ranges of ex- 
pected mass ratios and spin magnitudes, and for all spin 
orientations, to compute the probability of a given recoil 
speed. For example, for mass ratios between 1 and 10, 
and spin magnitudes of 0.9, the probability of exceeding 
1000 km s _1 is predicted to be ~ 10%. Studies of var- 
ious mo del- depend ent speed probattilities are given by 
ISchnittmanl (|2007l) and lBaker eToZI (|2008al ). 



2. Consequences of Black Hole Recoil 

Among the significant astrophysical consequences of 
gravitational-radiation recoil, growth rates of black holes 
can be affected. For example, the recoil of massive black 
holes imparting less than the escape velocity of their host 
galaxies tends to have regulatory effects on the amount 
of mass accreted ( Blecha and Loebl . |2008|) . Recoil ve- 
locities that exceed galactic escape velocities can im- 
pact massive black-hole growth in a differe nt way, by re- 
ducing the chances of subsequent mergers ( Sesana et all . 
120071: IVolonteril . 120071 ). This in turn may modestly re- 
duce the rate of coa lescence events observable by LISA 
(jSesana et ^.1 . 120071 ). 

The recoil might also have more directly observable 
consequences. Some quasars are thought to originate 
from the coalescence of massive black holes during galac- 
tic mergers. Quasars ejected from thei r host galaxies 
might therefore be expected. A study bv lBonning et all 
(l2007h found no evidence for such events, indicating 
that they are rare. However, recent o bservations of 
two r apidly moving extragalactic quasars (IShields et all . 
I2009D are strong candidates for such recoiled black holes 



(Komossa et al 



s candr 
., 2008) 



B. The Spin of the Final Black Hole 

The final spin of the merged remnant of a binary is also 
of astrophysical interest. It would be useful to know the 
probability of various spin magnitudes in constructing 
matched filtering templates and estimating gravitationa l 
signal detectability for LIGO and LISA ( Berti et all . 
l2007al) . The ability to predict a final spin given initial 
binary parameters also implies that observation of the 
spin of a black hole may help understand its origin. In 
particular, an understanding of the relationship between 
the final spin of a black hole and the orbital angular mo- 
mentum of its binary pre cursor may help explain the fo r- 
mation of X-shaped jets ( Barausse and Rezzollal . l2009h . 

As with the kick, it would be advantageous to con- 
struct a simple analytic formula with which to estimate 
the spin. Observations from numerical simulations sug- 
gest this might be possible. For example, in the case 
of nonspinning, equal-mass binaries undergoing circular 
inspiral, the final spin is insensitive to the initial separa- 
tion, being determined by the terminal dynamics of the 
merger . Such evolu t ion results in a universal final spi n of 
~ 0.7 dBaker et all . l2006al ICampanelli et all l2006allbl ldl: 
iPretoriusl l2005arK Additionally, for unequal-mass, non- 
spinning binaries, the final spi n scales roughly linear ly 
with symmetric mass ratio ([I} ([Gonzalez et all l2007bl) . 

Various approaches to constructing an ansatz for the 
final spin hav e been proposed. For alig ned-spin, equal- 
mass mergers. ICampanelli et all (|2006dl) produced a sim- 
ple formula quadratic in the (common) initial spin; this 
satisfied the "cosmic censorshi p" hypothe si s: af/ Mf < 
1. For nonspinning mergers, iBerti et al. I (l2007bh sim- 
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ply assumed a//Mf — ai) + br/ 2 , for fitting parame- 
ters a and b, which agreed with data from n o nspin - 
ning binaries reasonably well. lLousto et al\ ( 20101 ) 
arrived at a more generic ansatz, for initial black 
holes of arbitrary spin s , usin g the PN approximation. 
iBarausse and Rezzollal ( 20091 ) proposed a different ex- 
pression for generic binaries, using a set of assump- 
tions about the approximate conservation of the mag- 
nitudes and relative angles of certain angular momen- 
tum vectors. Several other formul a s have been suggested 
dBuonanno et all 120081: iRezzollal l2009t iRezzolla eiaH 



l2008allbl ; iTichv and Marronettii l2008j) . each in good 



agreement with some subset of available numerical data. 
However, broad agreement on an analytic formula for the 
final spin from gencrically processing binaries has yet to 
be achieved. 



C. Electromagnetic Counterparts of Black Hole Mergers 

We have seen above that black-hole mergers are 
expected to be "loud," since they produce strong 
gravitational-wave signals. But will they also be 
"bright?" That is, will there be an accompanying dis- 
play of photons, detectable by telescopes in any frequency 
range of the full electromagnetic spectrum? 



1. Astrophysical Considerations 

The answer to this question depends critically on the 
amount and distribution of gas and magnetic fields sur- 
rounding the merging binary. For stellar black-hole bi- 
naries, and intermediate-mass black-hole binaries in stel- 
lar clusters, any matter in an accretion disk would be 
consumed by one or both black holes and disappear rel- 
atively quickly, making an electromagnetic counterpart 
highly unlikely. 

Massive black-hole binaries formed from galaxy merg- 
ers, however, present a very different situation. In a gas- 
rich or "wet" merger, there is likely enough gas available 
to feed accretion disks around each black hole that even- 
tually evolve into a circumbinary disk; this provides a 
source of gas and magnetic fields that could generate de- 
tectable electromagnetic emission. Even in the case of 
gas-poor or "dry" mergers, the thin hot gas present is 
such (elliptical) galaxies might be sufficient to produce 
an electromagnetic signature. 

Such electromagnetic signatures would be valuable to 
astrophysics. Identification of the source on the sky 
would help confirm and characterize the merger, and 
probe accretion physics. A measurement of its red- 
shift using electromagnetic radiation, taken together with 
the determination of luminosit y distance using gravita - 
tional waves observed by LISA ^Lang and Hughesl . l2009l) . 
would provide an independent calibration of the distance 
scale and a precise probe of cosmolo gy including the na- 
ture of the mysterious dark energy ( Arun et ai , l2009bt 



Dalai et all 120061: iHolz arid Hughesl. 120051: [jonsson et all 
2007t IKocsis et all l2006ll2008ll2007t) ~Differences in ar- 
rival times of the electromagnetic and gravitational- wave 
signatures could also test fundamental principles such as 
the relative propagation speed of photons and gravita- 
tional waves. 

With the recent successes in numerical-relativity 
modeling of black-hole mergers, we find considerable 
interest in understanding possible electromagnetic 
signals from these events. Most work to date has 
focused on mechanisms that could produce emis- 
sion from a surrounding accretion disk, including 
signals induced by a recoi ling merged black hole 
encou n tering the disk (lArmitage and Nataraian , 
20021: iBode and Phinnevl l2007t IChang et all I200S; 



Megevand et all I2009T 



20051: 



O'Neill et al 



Dotti et all 120061: IHaiman et all [2009allbl: ; IKocsis et al 
20061 120081: IKocsis and Loebl. l2008t iLippai et all I200S ; 



.Milosavlievic and Phinnev 
20091 iPhinnevl 12007 



Schnittman and Krolikl . l2008l IShields and Bonning! 
2008[ ). Moreover, explorations of possible electromag 



netic signals that could arise in the dynamic spacetime 
near the binary during its last few orbits and merger are 
just now beginning. 



2. Simulations with Magnetic Fields or Gas near the Merging 
Holes 



iPalenzuela et al\ ( 2009f ) recently studied the effects of 
a merging black-hole binary on a surrounding magnetic 
field. The equal-mass, nonspinning holes start out with 
separation sa 6M on quasicircular orbits somewhat out- 
side the ISCO. The magnetic field is initially poloidal 
and assumed to be generated by currents in a distant 
circumbinary disk located at ~ 10 3 M. The electric field 
around the binary is initially zero. There is no matter 
near the black holes, as predicted when the circu mbinary 
disk is thin (jMilosavlievic and Phinnevl 120051) . They 
solved the coupled Einstein-Maxwell equations for a bi- 
nary evolving in presence of externally sourced electro- 
magnetic fields. As the system evolves, the inspiralling 
black holes stir up the fields. Figures I4U1 and |4"T1 show the 
electric and magnetic field lines at « 40M and « 20M 
before the merger, respectively. The magnetic fields are 
mostly aligned with the z axis, while the electric fields 
are twisted around the black holes. The binary dy- 
namics induce oscillations in the electromagnetic energy 
flux, with a period approximately half of the dominant 
(quadrupole) gravitational-wave signal. The energy in 
t he electromagnetic field is also enhanced gradually. 

IPalenzuela et all ( 20091 ) have certainly made an inter- 
esting start on this important problem. Can this sce- 
nario generate detectable electromagnetic emission? The 
answer awaits more realistic simulations, including spin- 
ning black holes. 

Other astrophysical scenarios allow both matter and 
magnetic fields in the vicinity of the binary close to 
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FIG. 40 Magnetic and electric fields lines around inspiralling 
equal-mass, nonspinning black holes at ~ 40M before the 
merger. from lPalenzuefa et ali l|2009f ). The electric field lines 
are twisted around the black holes, while the magnetic field 
lines are mainly aligned with the z axis. 




FIG. 41 Same as Fig. HOI except at a later time, « 20M 
before merger. 



merger. In this case, calculation of possible electromag- 
netic signatures requires solving the equations of gen- 
eral relativistic magnctohydrodynamics in a dynamically 
evolving spacetimc governed by the Einstein equations. 

Ivan Meter et al. I <I2010IJ) took a step towards solving 
this problem by mapping the flow of pressureless matter, 
modeled as non-interacting point particles, in the dynam- 
ical spacetime around the merging black holes. They 
started with a distribution of particles around the black- 
hole binary, and then evolved the binary using numerical 
relativity while tracing the motion of the particles along 
geodesies as the binary evolves. To estimate the energet- 
ics of the flow, they detected "collisions" by looking for 
particles within a small distance (typically, < 0.1M) of 
each other and then calculating the Lorentz factors 7 co u 
in the center-of-mass frame of each collision. 



They began with equal-mass inspiralling black-hole bi- 
naries ~ 5 orbits before merger, and considered both 
black holes to have either zero spin, or (a/M)\_2 = 0.8 
with spin vectors aligned with the binary orbital angu- 
lar momentum. A single black hole of mass M is also 
evolved as a control case. Approximately 75,000 geodesic 
particles are then initially distributed uniformly through- 
out a solid annulus centered on the binary and having 
inner radius 8M, outer radius 25M, and vertical full 
thickness lOAf. Particles within the inner radius are ex- 
cluded, to avoid transient signatures from particles ini- 
tially near the horizons. Note that these circumbinary 
disks are geometrically thick, and thus potentially have 
high enough inward radial speeds to keep up with the 
shrinking separation of the inspiralling binary, providing 
a source of gas near the black hol es during their merger 
( Milosavlievic and Phinnevl 120051 ) . 

Ivan Meter et al\ ( 2010b[ ) studied two initial velocity 
configurations for the particles. In the "orbital" con- 
figuration, the initial velocities are randomly distributed 
around a tangential velocity V c that would give a circular 
orbit in a Schwarzschild spacetime of mass M , resulting 
in a scale height of 5M. In contrast, the "isotropic" con- 
figuration is an extreme (that is, very hot) case in which 
the particles only have random velocities, with each com- 
ponent sampled from a Gaussian distribution of standard 
deviation V c j\fi. 

For the orbital initial velocities, the particles remain 
mostly in a rotational configuration; some particles do 
enter the region around the black holes, but are soon 
ejected by a gravitational slingshot. The collisions typi- 
cally occur at relatively low velocities with Lorentz fac- 
tors 7 C oii ^1.8. By the time of merger, there are essen- 
tially no particles near the black holes, and thus there 
is no energetic signature of the merger. The situation is 
quite different with the isotropic initial conditions, which 
provide a continual influx of particles throughout the evo- 
lution. Moreover, there is a clear signature of the merger 
visible in the maximum Lorentz factor, which takes on 
values 7coll,max ~ 2 during the inspiral and spikes up 
to 7coii.max ~ 3.5 just before merger for the nonspin- 
ning case. For rotating black holes, 7coli,max ~ 3 during 
the inspiral, spiking to 7 co li,max ~ 6 just before merger. 
In addition, gravitational torques from the merging bi- 
nary effectively stir the particle distribution, leading to 
high- velocity outflows with particles reaching 7 co u ~ 4 or 
more. 

Ivan Meter et al\ ( 2010bl ) point out that, in realistic 
astrophysical disks, viscosity would cause angular mo- 
mentum transport, bringing material inward towards the 
merging black holes. This could produce a scenario be- 
tween the two extremes the y studied, with a cl ear merger 
signal. More recent work bv lBode et all ll 2 01 Oft , using hy- 
drodynamical simulations of gas around black-hole bina- 
ries, provided strong support in this regard. Further con- 
firmation of this interesting possibility, along with other 
mechanisms involving general relativistic magnctohydro- 
dynamics, awaits more detailed simulations with gas and 
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magnetic fields, currently in development. 



IX. FRONTIERS AND FUTURE DIRECTIONS 

Five years ago, we did not know what obstacles might 
conspire to prevent the effective application of numerical- 
relativity simulations in deriving predictions of general 
relativity to address questions of astrophysical black-hole 
binaries. Up to that time, the numerical relativity re- 
searchers addressed their research primarily to questions 
of theoretical and computational physics. 

Then, quite suddenly, the last obstacles vanished. As 
reported, numerical relativity can now be applied to un- 
derstand general relativity's predictions of strong-field 
gravitational physics and to begin addressing questions 
of astrophysics and gravitational-wave data analysis. We 
expect future work in this field to be guided not by in- 
ternal problems in gravitational theory, but by questions 
of astrophysics, and other areas where strong-field grav- 
itational theory applies. These questions will motivate 
both continued, more detailed, investigation of some of 
the phenomena which have already been revealed, and 
the development of new capabilities to bring simulation 
results to bear on new questions. 



A. Gravitational-Wave Astronomy 

Though the work has begun in applying numerical- 
relativity results to gravitational-wave observations, 
much more work is still needed. It is now clear that 
the signal-to-noise (SNR) ratio from the last moments of 
the merger will dominate over the inspiral signal in many 
potential observations. Though numerical-relativity sim- 
ulations have shed considerable light on the basic fea- 
tures of the burst of gravitational radiation that com- 
pletes the predicted signals, applications in data analysis 
will require comprehensive quantitative knowledge of the 
merger signals generated over the full parameter space of 
mergers. 

While we discussed examples of mergers that sample 
the binary black-hole parameter space along what may be 
its principal full quantification of the signal space 

will require a long systematic program of investigation. 
Even with the aid of empirical models for encoding nu- 
merical results, many simulations must be conducted to 
qualify possible models. Initial models suitable for de- 
tection and (perhaps) parameter estimation of low-SNR 
signals, will require much further development for appli- 
cation to the high-SNR observations expected with fu- 
ture instrument upgrades, such as Advanced LIGO, and 
future instruments such as LISA. Considerably more sim- 
ulations of high quality will be required to achieve these 
goals. 

In addition to understanding the signals generated by 
mergers through the bulk of parameter space in greater 
depth, relativists must also extend the region of pa- 



rameter space covered by simulations. Typical cur- 
rent simulations arc limited in the range of parameters 
that can be practically covered. Currently popular ini- 
tial data models, particularly those assumin g conforma l 



flatness, are limited to (a /M) 1,2 ^ 0.93 ( Pain et all 



120021 ; lLovelace et aU l2008V Alternatives more suited to 
rapidly spinning Kerr black holes are less well developed 
in other ways. While astrophysical limits on spin mag- 
nitudes remain unclear, it is likely that larger spins will 
have to be treated with possibly novel initial data mod- 
els. Other, related challenges in the simulations may also 
exist. 

Current simulations are also restricted to comparable 
mass ratios. Simulations beyond q ~ 5 require some- 
what heroic investments of computational resources with 
present techniques. Astrophysical mergers may occur, 
however, over a broad range of mass ratios. In the limit 
of extremely large mass ratios, approximation techniques 
are possible that treat the motion of the smaller object 
as a perturbed geodesic in the spacetime of the larger 
black hole. Such analytic methods look promising in 
this limit, where numerical simulations are less practi- 
cal. However, there is a rather large middle ground, say 
10 < q < 100, where new techniques may be required. 
One problem is that the small spatial scale of the smaller 
black hole requires extremely small timesteps for stable 
explicit time integrations, even though there is nothing 
interesting happening on these tim e scales . Alternatives 
such as implicit schemes ( Lau et all . I2009T ) may open the 
door to a broader range of applications. 

Most of the simulations discussed also focus on circu- 
lar, or nearly circular orbits. While this is clearly an im- 
portant portion of the black-hole-binary population, sce- 
narios have been proposed in which one black hole cap- 
tures another from nearly parabolic initial encounters. 
Understanding the potential signals from such systems 
may require new techniques appropriate for long periods 
of slow, effectively Newtonian evolution, punctuated by 
brief periods of strong gravitational interaction. 



B. Other Astrophysics 

As discussed, the numerical discovery of strong 
gravitational- wave recoils in asymmetric mergers of spin- 
ning black holes has had significant impact in areas of as- 
trophysics beyond direct gravitational- wave observations. 
These results have spurred excitement in the black-hole 
astronomy community about other possible applications 
of numerical relativity. Since other areas of astronomy 
are based on (primarily) electromagnetic observations, 
these interactions naturally lead to applications involv- 
ing the interactions of black holes and electromagneti- 
cally visible matter. There are many relevant phenomena 
of interest: accretion disks around black holes, possibly 
disturbed by gravitational recoil, black-holc-neutron-star 
binaries, neutron-star-neutron-star binaries, jets from 
active galactic nuclei, quasars, and the mysterious ori- 
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FIG. 42 Numerical simulations are applied to study the crit- 
ical transition separating mergers fr om scattering events i n 
high- velocity black- hole encounters (jSperhake et all . l2009t ). 
The curves show the path of one black hole in each of three 
simulations begun with different initial conditions near the 
critical impact parameter. The trajectories track closely to- 
gether coming in from the right until the black holes encounter 
each other near the origin. Figure kindly provided by U. Sper- 
hake. 



gins of gamma ray bursts. 

A key direction for future work has been adding physics 
beyond the purely gravitational simulations we have re- 
viewed here. As discussed, efforts are currently underway 
to use and develop magnetohydrodynamics in a dynam- 
ical spacetime. As such simulations are applied with in- 
creasing realism to more astrophysical scenarios, refine- 
ment of the techniques will be required. Adequately re- 
solving shocks and turbulence interacting with magnetic 
fields, achieving robust accuracy in the presence of very 
strong magnetic fields, and enforcing the vanishing diver- 
gence of the magnetic fields, all in curved spacetime, are 
among the challenges that need to be met. 



C. Other Physics 

Astrophysical mergers typically involve nearly circular 
orbits, or in extreme scenarios, initially parabolic encoun- 
ters. Numerical simulations can be applied to other sorts 
of black-hole interactions as well. 

For high-velocity black-hole encounters, one expects 
a phase transition dividing the spacetimes in which the 
black holes merge from those in which they pass each 
other hypcrbolically, never to approach each other again: 
see Fig. 1421 Numerical simulation studies are begin- 
ning to elucidate c ritica l behavior near this phas e tran- 
sition (iPretoriusl. 120061: iPretorius and Khuranat l2007t 
ISperhake et aZ.1 . 12009^ . 

These fascinating phenomena at the extreme limit of 
strong gravitational dynamics have been discussed in 



the context of upcoming experimental measurements. 
In the trans-Planckian energy limit, some particle colli- 
sions may be describable by classical black-hole dynamic s 
([Banks and Fischlerl . ll999t lEardlev and Giddingsl . l2002t) . 
In highly speculative theories involving large extra 
spatial dimensions, it is conceivable that TeV scale 
physics may b e sufficien t to p r oduce such black- hole- 
like collisions (jGiddingsl . l200lfc iGiddings and Thomasl . 
l2002h . These possibilities motivate the first numerical- 
simul ation studies of ultra-relativistic black-hole coll i- 
sions (jShibata et all . 120081: ISperhake et all I2008U 120091) . 



D. Strong Gravity as Observational Science 

Einstein gave us general relativity nearly 100 years 
ago. Over the past century this theory, our standard 
model of gravitational physics, has passed experimen- 
tal and observational tests over ranging from laboratory- 
scale physics experiments and solar-system tests to obser- 
vations of c ompact astronomical objects and cosmology 
()Willl . 120061 These measurements have, so far, provided 
no need for a refinement of Einstein's theory of grav- 
ity [though some models to expl ain cosmic acceleration 
do in volve alternative gravity ( Silvestri and Trodden! . 

Numerical simulations are now revealing the detailed 
predictions of Einstein's theory for the final merger dy- 
namics of a black-hole binary, and its record in the emit- 
ted gravitational radiation. Gravitational-wave obser- 
vations of these events will expose such gravitational 
phenomena to measurements in a strong-gravity regime 
far beyond anything which has previously been tested. 
These include the first tests of higher-order terms in the 
PN expansion and, indeed, of the physical predictions 
which numerical relativity is just now revealing. 

The possibility that general relativity will find its limits 
in these observations motivates a better understanding of 
what wavefor ms might be predicted by alt ernative theo- 
ries of gravity ( Yunes and Pretoriusl [20091 ) . Alternative- 
gravity models of black-hole mergers would also likely 
require numerical simulation to derive waveform predic- 
tions, but little work has been d one in this direction so 
far [see one related example from lSalgado et al. I (120081) ]. 



If Einstein's theory proves correct in predicting the de- 
tailed physics revealed by numerical relativity, it would 
stand as a truly incredible achievement of scientific in- 
duction, divining the details of phenomena vastly re- 
moved from physical observations on which the theory 
was founded. If not, then these observations, together 
with a confident understanding of Einstein's predictions 
founded in numerical relativity simulations, may indeed 
lay the foundation for the next theory of gravity. 
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